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SUMMARY 


In this research project, an efficient multilevel adaptive method has been successfully applied to perform 
direct numerical simulation (DNS) of flow transition in 3-D channels and 3-D boundary layers with 2-D and 3-D 
isolated and distributed roughness elements resolved by a curvilinear coordinate system. The approach consists of 
a fourth-order finite difference technique on stretched and staggered grids, a fully- implicit time marching scheme, 
a semi-coarsening multigrid method based on line distributive relaxation, and an improved outflow boundary- 
condition treatment, which needs only a very short buffer domain to damp all order-one wave reflections. This 
approach makes the multigrid DNS code very accurate and efficient. It allows not only spatial DNS for the 3-D 
channel and flat plate at low computational cost, but also spatial DNS for transition in the 3-D boundary layer 
with 3-D single and multiple roughness elements. These calculations would have extremely high computational 
costs with conventional methods. Numerical results show good agreement with linear stability theory, secondary 
instability theory, and a number of laboratory experiments. The contribution of isolated and distributed roughness 
to transition is also analyzed. 


1 Introduction 

The main driving force for the study of flow transition is the understanding, prediction, and control of transition 
to turbulence. The transition process from laminar to turbulent flow in a shear layer is a basic scientific problem 
in modern fluid mechanics and has been the subject of study for over a century, but it remains a challenging 
and important area for research (Reed & Saric, 1989). The control of flow transition to turbulent state is a 
very significant topic. The drag of an aircraft at cruise flight conditions is about 60% friction drag for present- 
day transport aircraft with turbulent boundary layers on their wetted surface. The principal drag reduction 
opportunities lie in stabilizing the laminar boundary layer as much as possible (Reshotko, 1988). However, the 
development of flow transition control is mainly based on our knowledge of flow transition. Therefore, developing 
an understanding of shear flow at high Reynolds numbers has been a central problem in the theory of fluid motion 
(Bayly, Orszag, fe Herbert, 1988). 

Classical linear stability theory, which was established independently by Orr (1907a, 1907b) and Sommerfeld 
(1908), is based on a parallel base flow assumption on the linearized perturbation form of the Navier-Stokes 
equations, referred to as the Orr-Sommerfeld equation. Tollmien (1931) and Schlichting (1932) were able to 
solve the Orr-Sommerfeld equation about 20 years later. Linear stability theory validation was first confirmed by 
Schubauer and Skramstad (1948), who used a vibrating ribbon to impress a disturbance into the boundary layer 
and hot wires to take measurements. Linear stability theory is now widely accepted as a tool to predict the initial 
growth of disturbances, more commonly referred to as Tollmien- Schlichting (T-S) wave, which can be obtained 
by solving the Orr-Sommerfeld equation (Joslin et al., 1992). 

When the amplitude of a T-S wave reaches a certain finite value ( about 1% of the base flow maximum 
velocity in root-mean-square (rms) values ), nonlinear effects are felt and the instability waves no longer obey 
the linear theory. At this stage, the three-dimensional disturbances are energized and streamwise vortices form 
and intensify. Later, instantaneous streamwise velocity profiles become inflectional, resulting in high shear layers, 
which is accompanied by the appearance of A-like structures. Eventually the high shear layers develop kinks, 
leading to the formation of hairpin eddies followed by the breakdown of the flow. Finally, turbulent spots form, 
leading to turbulence. 

Nonlinear instability, commonly referred to as secondary instability, is associated with the formation of three- 
dimensional structures. An aligned three-dimensional structure was first observed by Klebanoff, Tidstrom, & 
Sargent (1962), now referred to as fundamental or K-type after Klebanoff. The staggered patterns of three- 
dimensional structure characterized by subharmonic signals were found by Kachanov, Kozlov, & Levchenko 
(1978). These experiments showed that the subharmonic of the fundamental wave (a necessary feature of the 
staggered patterns) was excited in the boundary layer. At small amplitudes, it produced the resonant wave 



2 


Liu, Liu, McCormick, Final Report 


interactions predicted by Craik (1971). This type of instability is referred to as C-type after Craik. At large 
amplitudes, Craik’s mechanism becomes less important and instability is characterized by the dominance of off- 
resonance modes. This type of instability is referred to as H-type after Herbert (1983a). Benney and Lin (1960) 
imposed a two-dimensional T-S wave and a pair of oblique three-dimensional wave at the inflow boundary to study 
the interactions between them. Later, Herbert (1983a, 1983b) derived a theory for an experimentally observed 
three-dimensional parametric instability. Remarkable agreement between prediction from this new theory and 
experimental results is obtained, and the secondary instability theory is now generally accepted as a tool to 
predict later development of transition in boundary layers. 

The major engineering tool for transition prediction is the e N - method. Smith and Gamberoni (1956) corre- 
lated transition data with the amplitude ratio of the most amplified eigenmodes between the onset of instability 
and the transition location and found that transition occurred when the amplitude ratio reached a value of 
e 9 « 8100. However, the method is semi- empirical and thus requires some foreknowledge of the flow undergoing 
transition. 

Although great progress has been achieved in this century, the mechanism of transition remains a mystery. 
Flow transition is a very complicated, three-dimensional, time-dependent, multistage physical phenomena, and 
is affected by many factors, including disturbance environment, pressure gradient, surface roughness, three- 
dimensionality, curvature, Mach number, and many other aspects that influence boundary layer development. 
For example, the manner of the effect of surface roughness on flow transition is still largely unknown. This is one 
of the major interests of the current work. 

The appearance of modern digital computers led to the development of numerical methods applied to studying 
instability and transition from laminar to turbulent flow. The most successful approach is by direct numerical 
simulation (DNS). Due to insufficient computer resources, the majority of the first investigations were temporal. 
Among them are those by Orszag & Kells (1980), Wray Sc Hussaini (1984), Biringen (1984), Kleiser Sc Laurien 
(1985), Zang Sc Hussaini (1985), and Zang, Krist, Erlebacher, Sc Hussaini (1987). Temporal simulations follow 
the time evolution of a single wavelength of the disturbance as it convects with the phase speed of the wave. 
Available mesh resolution can then be focused on resolving a single wavelength. This enabled simulation into 
the later stages of transition and revealed good agreement with the experiments. However, the assumption of 
periodicity along the streamwise direction does not allow direct comparison of temporal numerical simulations 
with spatially evolving experimental results. Moreover, streamwise growth of the base flow (e.g., in boundary 
layers) is not allowed in the temporal approach. Hence, boundary layer simulations are restricted to the parallel 
flow assumption as well. This approach basically lacks a physically realistic representation. 

On the other hand, the spatial DNS approach provides the needed quantitative information about transition. 
But with spatial DNS, obstacles have existed that have prevented fully carrying out such a study. Among these 
are the realistic specification of inflow and outflow boundary conditions and the high demands on computational 
resources. Even with today’s supercomputers, current resources are insufficient to allow conventional methods to 
perform full simulation of transition to turbulence in a boundary layer in a spatial setting. 

A number of authors have made some progress in spatial DNS, including Fasel (1976), Fasel Sc Bestek (1980), 
Fasel Sc Konzelmann (1990), Spalart (1989), Danabasoglu, Biringen, Sc Streett (1991), and Joslin, Streett, Sc 
Chang (1992). The results provided by spatial DNS show qualitative agreement with linear stability, secondary 
instability theory, and some available experiments. Kleiser Sc Zang (1991) gave a more complete list of DNS works 
in their review paper. It is worth mentioning that significant achievements in spatial DNS have been made by Liu, 
Liu, Sc McCormick (1991, 1992, 1993). Supported by the NASA Langley Research Center since 1990, they have 
developed a very efficient and accurate multigrid approach for DNS. The flow transition process, including linear 
evolving secondary instability and breakdown in 3-D channels and 3-D flat plates, was fully simulated, showing 
good agreement with linear stability theory, secondary instability theory, and experiments by Nishioka, Asai, Sc 
Iida (1981) and Saric Sc Reed (1987). Since the approach is computationally very efficient, it enabled simulation 
of more complicated transition processes such as with 3-D distributed roughness elements in 3-D boundary layers. 
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Instability of boundary layers with surface roughness is an important topic because of its practical applications 
to describe, model, and control transition for efficient airplane design. Several experimental studies have been con- 
ducted in boundary layers with two- and three- dimensional isolated or distributed roughness elements. Klebanoff 
k Tidstrom (1972) used a cylindrical trip wire as an isolated roughness element and showed that roughness itself 
does not introduce disturbances, but instead causes early transition via the destabilising influence of the flow in 
the recovery zone. These disturbances rapidly amplify in the streamwise direction and cause the transition point 
to move more upstream. Boiko, Dovgal, Kozlov, k Shcherbakov (1990) and Dovgal k Kozlov (1990) used a rect- 
angular two-dimensional isolated roughness element, and excited the boundary layer at several frequencies by a 
vibrating ribbon. These investigation revealed the growth of instability with frequencies in the inviscid instability 
range, showing that the development of disturbances in the recirculation zone is governed by the instability of the 
mixing (shear) layer at its edge. A more realistic case that has distributed roughness elements in a 3-D channel 
has been experimentally studied by Tadjfar, Reshotko, Dybbs, and Edwards (1993) most recently. They found 
that the oscillation is amplified rapidly as it moves downstream, but there is no evidence of Tollmien-Schlichting 
instability waves in the flow. 

On the other hand, existing numerical studies of surface roughness are very limited. The only computational 
studies concerning the spatial stability of a separated flow were conducted by Fasel, Bestek, k Schefenacker (1977) 
and Danabasoglu, Biringen, k Streett (1993). Fasel et al. (1977) investigated a two-dimensional numerical sim- 
ulation of a boundary layer with a backward-facing step. Danabasoglu et al. (1993) investigated an isolated, 
two-dimensional roughness element on the spatial development of instability waves in 2-D boundary layers. Ap- 
parently, all existing numerical studies on surface roughness are limited to a single 2-D roughness elements in a 

2- D boundary layers. After successfully simulating the flow transition from linear instability to breakdown in 3-D 
channels and 3-D boundary layers in a very efficient way (usually it takes less than 6 CRAY hours for each case), 
Liu k Liu continued their study and have successfully developed a fast multigrid solver, which can now solve 
the 3-D, time-dependent, incompressible equations in general coordinates. This solver has successfully simulated 
transitional flow past multiple 2-D and 3-D roughness elements with general roughness shape and distribution in 

3- D channels and boundary layers. 

The new approach used in this work includes: 

• General coordinates that are able to work for more general geometry 

• Multilevel stretched and staggered grids that are generally accepted to have higher accuracy than that of 
non-staggered grids 

• Fully implicit time-marching (implicit for both convection and diffusion) that has much better stability than 
that of explicit methods 

• Fourth-order discretization in all three directions that has high accuracy and can reduce artificial viscosity 
and phase error 

• Line distributive relaxation that is very efficient smoother for Navier-Stokes equations 

• Semi-coarsening multigrid for solving the large scale nonlinear algebraic systems arising in the implicit 
approach at each time step based on line distributive relaxation to make fully implicit schemes comparable 
to explicit schemes at each time step in CPU cost 

• New out-flow boundary treatment that efficiently eliminates all wave reflections from the outflow boundary 
and reduces the size of buffer domain to as short as less than one T-S wavelength for both small and large 
magnitude of inflow disturbance, which greatly enhances the efficiency of spatial DNS 


Since the new approach developed in this work has high accuracy, stability, and efficiency, not only has 
spatially linear growth and secondary instability in 3-D channels and flat plates been successfully simulated in an 
efficient way, but also new numerical simulation regimes have been reached: successful simulation of transition 
over multiple 2-D and 3-D roughness elements, with general shape and general distribution in 2-D and 3-D 
channels and flat plates. 
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2 Governing Equations in Physical Space 


2.1 Original Form 

Low speed flows are governed by the incompressible Navier-Stokes equations, in which momentum and mass 
conservation are represented by three momentum equations and a continuity equation, respectively. Considering 
only Cartesian coordinates, the governing equations can be expressed as follows: 
x-momentum equation: 


y - momentum equation: 


z-momentum equation: 


continuity equation: 


Du* 

dP* 

P Dt* 

1 

Qsl 

H 

* 

Dv* 

dP* 

P Dt* 

dy* 

D- w* 

dP* 




u or A . * . 


Dt* 


Dp .du* dv* dw * . 

♦ + ^ + a7 + ^ ) = 0 - 


Dt 


( 1 ) 

( 2 ) 

(3) 

(4) 


where, u*, u*, in* are dimensional streamwise, normal, and spanwise velocity components, respectively, P* is the 
dimensional pressure, p represents the density of the fluid, /*, / y , and f x are body forces in the three directions, p 
is the dynamic viscosity coefficient, A is the three-dimensional Laplacian, and denotes the material derivative. 
In Cartesian coordinates, we have 


A* = 


d 2 


dx 


^ + 


d 2 

d y * 2 


4 


d 2 

dz* 2 


D d . d . d t . d 

— -+■ u 4- v — “I - in — — . 

DV dV dx* dy* dz* 


It is assumed that no density variation exists in the investigation domain, thus, 
and the continuity equation becomes 


Dt* 


du * dv * dw * 

dx* dy* + dz * 


0 . 


(5) 

( 6 ) 

(7) 

( 8 ) 


Furthermore, assuming the effect of the body force / = (/*, / y , /,) is negligible, the momentum equations can be 
written as 


du* 

dt* 


t du* 

dx* 


, du* 

dy* 

t dv* 

dy* 


. du* 1 dP * 

u* j 

dz* p dx* 


dv* . dv * .dv* .dv* 1 dP* 

dt* dx* dy* dz* p dy 


dw* . dw * . dw 

-I- u* — — + v ^r— '+ ™ 


dt* 


dx* 


dy* 


. dw* 1 dP* 
dz* p dz* 


f d 2 u* 


d 2 u* 


d 2 u* 


(9) 


4- 

dy * 1 

4- 

d. i* 1 ’ 

= 0, 

,d 2 v* 


d 2 v* 


d 2 v* 

= 0, 

(10) 


■ 4- 

dy* 1 

4* 

dz* 1 ’ 

, d 2 W* 


d 2 w* 


d 2 w* 


(11) 

ax- 

4- 

dy * 1 

4* 

dz * 1 * 

= 0. 


Here, v — is the kinematic viscosity coefficient. 
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2.2 Nondimensionalized Forms for Planar Channel-Type Flow 

Equations (8)-( 11) can be nondimensionalized by dividing them by some reference constants appropriate to 
the flow, for instance, reference velocity and reference length. For planar channel-type flow, we use the channel 
half height h as the reference length and the center line velocity Uo of the steady flow as the reference velocity. 

Define the nondimensional variable as: 


and the Reynolds number 



u* 

v" 

U — 

Uo’ 

V ~Uo’ 


X* 

y* 

X — 

T’ 

y= -h' 


p* 

, CU 0 

p = 

W’ 

h 



Re 


U 0 h 

j 

v 


we finally obtain the following nondimensionalized governing equations: 


du du du du dP 1 .d 2 u d 2 u d 2 u . 

_i_ ii — 4 - v — 4 - w 1 — ( 1 j — ) — 0 , 

dy dz ^ dx Re K dx 2 ^ dy 2 ^ »- 2 } 


dt 


dx 


dv d v d v dv dP 

1- u h v h w (- 

dt dx dy dz dy 

dw dw dw dw dP 

L u L V L w ^ _ — 

dt dx dy dz dz 


dz 2 J 

d 2 v , 


1 d 2 v d 2 v v 

Re dx 2 + dy 2 dz 2 1 

1 d 2 w d 2 w d 2 w 

Re dx 2 dy 2 dz 2 ’ 

du dv dw 
dx ^ dy ^ dz 


( 12 ) 

(13) 

(14) 

(15) 

(16) 
(17) 


It is recommended to use the conservative form of the momentum equation in the numerical simulation, which 
can be derived from (14)-(17): 


du duu duv duw dP 

dt ^ dx ^ dy ^ dz dx 

dv duv dw dvw dP 

dt dx + dy + dz ^ dy 


1 d 2 U (Pu d^u _ 
Re dx 2 + dy 2 + dz 2 


J_ d 2 ^ d 2 ^ d 2 ^ 

J2e dx 2 ^ dy 2 dz 2 1 

dw duw dvw dww dP 1 d 2 w d 2 w d 2 w 

H o I o t ^ I" d 71 a. .2 q z 2 / — * 


dt 


dx 


dy 


dz 


dz Re K dx 2 dy 2 


(18) 

(19) 

( 20 ) 


The major purpose of this work is to investigate the behavior of disturbances, so the perturbation form of 
the equations plays a more important role. The perturbation equations are obtained by decomposing the total 
flow into steady base flow and a fluctuating unsteady perturbation. Letting the subscript 0 denote the base flow 
variables and prime the fluctuating variables, we have 


«( z > y, 2 , 0 = y< z ) + u '( x . y> z > *)> 

v(x, y , z, t) = u 0 (z, y, z) + d'(x, y, z, t), 
tu(z, y,z,t) = w 0 (x, y, z) + w’(x, y, z, t), 

P(x, y, z , t) = P 0 (x, y, z) + P'(x, y, z, t ). (21) 


Substituting (21) into (17)-(20), and noting that the base flow itself also satisfies the Navier-Stokes equations, we 
have 
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sc-momentum equation: 


du l d{2uQv! 4 u'u') d(uov* 4 ti'vo 4 uV) 

+ dx + dy 

d{uow' 4 v!wq 4 tiV) dP' 1 d 2 u 1 d 2 u ' d 2 v! 

d~ z + ~dx ~ fc^-fc 2 + ~dy 2 + dz 2 } 


= 0 , 


( 22 ) 


y-momentum equation: 


dv' d(u'v 0 + uo^' + uV) d(2v 0 v' + v'v') 
~dt + dx + dy 

d{v 0 w' + v'wq + v'w') dP' 1 d 3 v' d 3 v' d 3 v' 
dz + dy Re dx 3 + dy 3 + dz 3 ' 


0 , 


(23) 


2 -momentum equation: 


dw ' ( d(u 0 u/ 4 u l wo 4 u'w*) d(v'w 0 4 v 0 w* 4 v'w*) 
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\3 f d 2 w ' d 2 w ' 

*■» o n 0 / 


d(2w'wo 4 w'w') dP 1 
dz dz 


Re ^ dx 2 T dy 2 r dz 2 


0. 


(24) 


continuity equation: 


du l dv' dw ' 
dx ^ dy ^ dz 


0. 


(25) 


For simplicity, we will omit all primes henceforth, so that all variables without any sub/superscript will denote 
the fluctuating part of the flow in the perturbation equations. In planar channel flow, for example, since both v 0 
and w 0 vanish and u 0 = u 0 (y), *he a ^ ove momentum equations can be simplified to 


dt 


dv 

dt 


4 


dt 


dx 

>(u 0 v 4 
dx 
w 4 uv 
dx 


dy dz 

dw dv w 

^ dy dz 

dvw dww 
dy + dz 


dP 

1 d 3 u 

d 3 U 


d 3 u. 


(26) 


~ ~Re'd^ 3 

+ dy 3 

4 
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= o, 

dP 

1 ,d 3 v 

d 3 v 


d 3 v . 

= 0, 

(27) 

+ 5y 

~ Re^dx 3 

+ dy 3 

4 


dP 

1 d 3 w 

8 3 w 



= 0. 

(28) 

h 37 “ 

' AT 

^ dy 3 

4 

dz 3 ^ 


2.3 Nondimensionalized Forms for Flat-Plate-Type Flow 

There are several methods to nondimensionalize the Navier-Stokes equations for flat-plate-type flow. Since 
the choice of reference length for this kind of flow seems to be somewhat arbitrary, we can for instance choose the 
distance from the leading edge to the inflow boundary, L 0} or the thickness of the boundary layer at inflow, 6q } 
or the thickness of the displacement at inflow, 6 q, or something else, as the reference length. In this paper, we 
choose Sq as the reference length, and the free stream velocity Uoo as the reference velocity. Under this setting, 
the Reynolds number is defined as 


Re' 0 = (29) 

V 

The general form of the governing equations is the same as (14)-(17) if we substitute Re by ReJ. For the 
perturbation flow, if the flat plate is used, the base flow is two-dimensional, and tu 0 vanished from equations 
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(22)-(24). Thus, the resulting governing equations 


are 


du d(uu -f 2uoii) d(uv -f uqv + uuo) 
dt dx dy 


diuw + uquj ) dP 1 ,d 2 u d 2 u d 2 u . 

" 8* Ret K ax? dy 2 dz 2) ’ 


dz 


dx 

dP 


dv d(uv + uou -j- u-uo) d(vv + 2vqv) d(vw + vow) 
dt ^ dx + dy + dz + dy 

dw 


Reydx 2 

1 d 2 v d 2 v d 2 v 

Re$dz 2 dy 2 + dz 2 
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dt 


^ d(uw + uow) ^ d(vw + votu) dww dP 1 ,d 2 w d 2 w d 2 w 


dx 


dy 


dz 


l/l A - 1/ vu 
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Re'^dx 2 
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) = 0 , 


£?u dv dw 

si + "57 + si - °- 


(30) 

(31) 

(32) 

(33) 


2.4 Initial and Boundary Conditions for Planar-Channel Type Flow 

For a planar channel, Poiseuille flow is considered as the base flow: 

U{y) = i-(»- i) J . ye [o,2]. (34) 

The initial conditions for the perturbation flow for temporal simulation are 

(u,v,id) = ( 2 dReal{(<l> u ,<f> v , 0) 3rf e* (<, * ) } 

+e 3 d + -Rea/{(^ UI Khd+e i{aX+t>,) } 

+^3 d- Real{{4>u,<t>v,K) 3 d-e i(ax ~ pl) }, 

P = 0, (35) 

where <f> Ui <f> v and <f> w are eigenfunctions for u , v , and u; obtained from Orr-Sommerfeld equation, 634, €34+ and 
€3*- are amplitudes of the 2-D and 3-D perturbation waves, oci and ot 2 are streamwise wave numbers for 2-D and 
3-D perturbation, and (3 is the spanwise wavenumber. For spatial simulation, the initial conditions are 

(u, v , w, P) = 0. (36) 

On solid walls, no-slip conditions are assumed: 

(u,v l ttf)| aa {i = 0. (37) 

No boundary condition for pressure is needed since we use a staggered grid for the numerical simulation. 

For temporal simulation, periodic boundary conditions are imposed in both the streamwise and spanwise 
directions: 

V(x + L x ,y,z) = V(x,y,z), 

P(x + L x ,y,z) = P(x, y, z), 

V(x,y,z + L t ) = V(x ,y,z), 

P(x,y,z + L,) = P(x,y,z). (38) 

For spatial simulations, the inflow boundary conditions are obtained from the 2-D T-S wave and a 3-D oblique 
wave pair. Because of the specialty of staggered grid, the condition can be written as: 

= €2dReai{4>u2de~ tu,t } + ^3 d+Real{4>u3d+e x ^ x - ut ^} 

+ e 3d _Real{4> u3 d-e^- p, - u ^}, 

= € 2 d Real{4> V 2de' l '~ a ^'~ ut) } + ( 3d+ Real{<i> v3d+ e^~ a ^ L+p, ~ ut ' 1 } 

+ € 3d _Real{<f> v 

= (3d+ Real{<l> w3d+ e'l- a ^ +p *-^} + (39) 


u(o,y, 

z,t) 

Ax 

z,t) 


Ax 


w ( — 7T'y> 

z,t) 
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The outflow boundary conditions are much more complicated, and will be discussed in detail in later sections. 

Periodicity conditions are also used in the spanwise direction for 3-D temporal and spatial simulations because 
the flow in this direction is not confined, and the laboratory experiments also support the spanwise periodic 
structures. 


2,5 Initial and Boundary Conditions for Flat-Plate-Type Flow 

For flows over a flat plate, the Blasius similarity solution is applied as the base flow. For rough plates, the 
Blasius solution will be used as the initial guess of the base flow. Details will be given in later chapters. 


Only spatial simulation is conducted for smooth or rough plates in this study. Thus, similar to the last section, 
Benny-Lin (1960) type disturbances are imposed at the inflow boundary: 


u(0,y,z,t) 





u>( — 


Ai 

T* 


y,M) 


= e 3d Real{<f> u3d e *"*} 

+ e 3d+ iZea/{^ + e’ (/, ‘- wt) } 

+ e 3d+ iiea/{^3a + e‘ ( - a ^ +<,I -“‘ ) } 

+ e 3d _Real{^> v3d -e^- a ^-^- u % 

= e 3d+ Real{<f> w3d +e i (~ a ^' + P‘- wt ^} 


(40) 


where w is the real frequency of the disturbance, /? is a real number that represents the spanwise wavenumber, 
and a = a/* + iaj is the streamwise complex wavenumber obtained from linear stability theory. 

A no slip boundary condition is applied at the solid wall. According to the linear stability theory, the 
disturbances vanish at infinity, so, we obtain the boundary conditions at far field given by 


u(x,y-+ oo,z,t) - 0, 
v(x,y-> oo ,z,t) = 0, 

w(x, y — > oo, z, t) = 0. (41) 

Also, no pressure condition is needed at the inflow, solid wall, and far field since a staggered grid is used. 


The outflow boundary condition will be discussed later. 


3 Transformation of Governing Equations 

3.1 Vector Form of the Navier-Stokes Equations 

The dimensionless 3-D incompressible Navier-Stokes equations in Cartesian coordinates are given in (17)-(20). 
It is convenient to rewrite them to a general vector form. First, we introduce the following vectors: 



1 


u 


V 


w 


u 


u 2 + P 

_ 

vu 

- 

wu 

u = 

V 

, E = 

uv 

, F = 

V 2 + P 

, G = 

wv 


w 


uw 


wv 


W 2 +P 


Liu,Liu } McCormick, Final Report 


9 


Then equations (l7)-(20) can be rewritten in a simple form: 

dV_ dE_ 5F dG _ J_ - 

dt ^ dx dy dz Re * 

where 

d 2 d 2 d 2 

A = 1 1- 

dx} dy 2 dz 2 

is the Laplacian operator in the physical (z,y, z) space. 


(43) 


3.2 General Transformation 

In order to treat arbitrarily shaped boundaries, it is often convenient to employ general curvilinear coordinates 
l, *?> C> which are related to Cartesian coordinates by a transformation of the general form 


* = C), 

y = y{l,v,0, 

or the corresponding inverse transformation 

l = ((*,y,z), 

r) = r)(x,y,z), 

C = 

Using the chain rule, the partial derivatives become 

a_ a a_ a 

dx~^ z di +T>x d^ +< ' x d(’ 

a a a a 

a^- iv dC +r,v a^ + Cv ac' 
a a a_ a 

8^ ~^d( +11l dr] +<:i d(' 


(44) 


(45) 


(46) 


If the inverse analytical transformation (45) is easy to obtain, the metric £ Xi rj X) £ x , 7] y , ( y , f 2 , rj ti ( 2 can be 

obtained directly. Otherwise, the following procedure is recommended. Rewriting the differential expressions in 
matrix form, we have 


and, similarly, 


dl ' 


1 

V? 


dx 

dr] 

= 

>7* ‘Hy Vi 


dy 

d( . 


. Cv Cx . 


dz 


dx 


X( X, x ( 


' ' 

dy 

— 

y( Vr, y< 


dr] 

dz 


. z ( Z V z < . 


. dC . 


Defining the Jacobian of the transformation J as 


dj^V, C) 
d{x,y,z) 


lx ly lx 

Vx Wy Vz J 

C* Cy C z 


(47) 


(48) 


( 49 ) 
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from (47) and (48) we obtain 


■ t. 


6 ‘ 



x v 

x < 

-1 


-{x„z ( - x^z,) 

- x (Vr, 

Vz 

% 


= 

V( 

Vv 

y< 

= j 

-( vt z < - y< z () 

X(Z( - x < z f 

~{ x iV( - X (V() 

. <* 

c. 

Cz . 


. H 

z v 

z < . 


yt 3 ^- V‘taZ( 

-(x f z, - x,z f ) 

x (Vv ~ x vV( . 


and, thus, 

tz = •%!*<- »<**). 

iy = — J(x v Z( — X(Z, ,), 

(m = J{ x r\VC ~ x (yt))> 

Vx = ~J(V(Z( — V( x (), 

Vy = •?{*(*( 

T) t = -j{x ( y ( - x ( y ( ), 

C* = J(y(*v - Vv z ()> 

Cy = ~J{ x t z V — x v z t)> 

Cm= J{ x <yv- X iy()- ( 50 ) 

After the proper transformation, a domain with general geometry can be transferred to a cubic domain in the 
computational (£, 77 , C) space. 


3,3 General Curvilinear Coordinates System 

Based on the above transformation, we now define the new vectors in the computational space: 


Vi = 


u 

7 ’ 


E 1 — + E'iy -f G^,), 

Fi ~ + Frjy + Grj s ) } 


G 


(43) then becomes 


l = y(fiC. + *C, + OG)- 


atfx dE 1 dF x dG x _ 1 - 

_ ar + 'ar + a^ + ~a< _ - ffe Al ^’ 


(51) 


(52) 


where Ai is the physical Laplacian operator transferred to the computational ((,r),Q space: 

„ d d , d . d d , d . 

Ai “ + ^ a^ + ^ d( +1]x di ] +< * x d(^ 

a a . a . . a a . a 

+( ^ a* + + c 'd( ){(y dt + ^ 

„ a a , a . a a . a , 

+(6 ^ + ^ + & 3^* + ^ a»? + ^ ac 


(53) 


Though (52) looks very simple, it is not convenient to discretize. To obtain equations that are more easily 
discretized for numerical simulation, we define the new variables 17, V , and W in the computational (£, 77 , C) 
space: 
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v = ji^Vx + vT} y + writ), 
W = Cx + Ky + WC), 


(55) 

(56) 


and consider both variable triples u, v, w and U, V, W as unknowns in the computational space. Our objective 
is to develop relatively simple expressions for the governing equations in general coordinates that are easy to 
discretize. Using the above transformation, we rewrite the vectors U\, E\, F, and G\ in the following form: 


Ul = J 


Fi 


1 

u 

11 

w 


V 

uv+?y± 

vV + 
wV+ 


E i 


i Gi 


r u 

uU + ^- 
vU + ^ 

L j 

w 

uW+ ^ 
vW + TT* 

L wW + fy J 


(57) 


The governing equations then become 


du . dUu dVu dWu 

ai + J{ ~ar + l^ + 


d( ^ + ^ X d( +r>x dri +< '*d(^ P " AlU 


dv 

at + J{ ~at 


dUv dV v dW v 


+ 


dr) 


+ 


) + (€. 


d 


d 


d( ' ' K " y + Vy dr, +< > y d^ P 


dw , dUw dVw dWw d d d . 

It + j (tt + + it’ + <f - Sf + ”■ T v + c - a< )p 


Re 

1 

Re 

1 

We 


Aiu 


0, 

0, 


Aiu; — 0, 


dU dV dW 

~at + + Ik 


= 0, 


(58) 

(59) 

(60) 
(61) 


The system (54)-(56) and (58)-(61) has 7 unknowns, u, v, w, U, V, W and P. Similarly, we can obtain the pertur- 
bation form of the governing equation as follows: 


du Jf d[u{u + Uo)+uoU] d[u(V + Vo) + Uo^] 
dt + J[ di + dr, 

a[u(vy-t-Wo) + u 0 w], ,d a . a , n 1 A 

— df = °’ 

dv Tf d[v{U + U 0 ) + v 0 U] 9[u(V + Vb) + v 0 V] 

at + ( at + a v 


B[v(w + Wo) + voW\ a_ a_ ± )p _± A 

d( ) + ^ y d( +Vy dr, + < ‘ y dC ) Re 1 

dw T/ d[w(U + U 0 )+ w 0 U] d[w(V -(- Vb) + w 0 V] 
at + J{ at + dr, 

a [u{ w + w 0 ) + woW] + ^± + n ± + (t ± )P _ J_ AlU; 
at ’ at dr, K at 1 Re 

au av aw 

~at + ~&n + ~at 


o, 


o, 


= o. 


(62) 

(63) 

(64) 

(65) 


where Uq, Vq and Wo are base flow solutions defined in the same way as (54)-(56). Combined with (54)-(56), this 
system also has 7 equations and 7 unknowns for the perturbations. From (53), we define the operator Aj in the 
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computational space as 


A x 


(ti +^+ 0^2 + +*$+ + (c« + c ; + <?>£ 

+ 2 (^ x 17 * -f CyWy + C*Vm) 2(£*C* + £*C*) 

+ 2(7fc £* + 77y^y + 7] X £*) Q^Q£ ^ (^** + Cyy "♦" ^X*)^ 

+ (*/** + *}yy + *7xx)^ H (C** + Cyy + Cm *)^T* 

^77 dC 


( 66 ) 


Because of the high accuracy requirement in flow transition simulation, an analytical transformation is recom- 
mended, so, there will be no extra artificial errors in the system (54)-(56) and (58)-(61) for the base flow, or 
(54)-(56) and (62)-(65) for the perturbation flow introduced by the transformation. Note that no orthogonal 
assumption for the curvilinear coordinates is made, so the above governing equations can in principle handle any 
kind of curvilinear coordinates, provided the transformation exists and is one to one. The new emphasis here 
is that the system has 7 unknowns (it, v, w, U, V ) W and P) instead of the usual 4 unknowns (17, V, W and P ) 
suggested by others. 


4 Numerical Methods for Two-Dimensional Problems 

The Navier-Stokes equations and their transformations given in section 2 and section 3 are nonlinear, time- 
dependent and three-dimensional. Except for a few cases, analytical solutions are impossible to obtain, and the 
numerical methods must be used. For numerical simulation of transition, there are several requirements we must 
meet: 

• phase-accurate discretization of the convective terms 

• high resolution for regions close to the solid wall 

• stability of the discrete system 

• proper outflow boundary conditions that remain non-reflective even in the presence of nonlinear wave 
interaction 

In addition to the above basic requirements, efficiency, memory storage, and vectorization of the computational 
program must also be considered. 

In this section, we first present an accurate method that meets these requirement for the two-dimensional 
problems. The procedure uses fourth-order central finite differences in space and second-order backward Euler in 
time. The numerical procedure is conducted on a staggered grid, so the resulting schemes are quite complicated 
to describe. At the outflow boundary, the buffer domain (cf. Streett fc Macaraeg, 1989) method is applied. 


4.1 Discretization on a Uniform Staggered Grid 

A uniform staggered grid is first employed to discretize the computational domain in the case of flow past a 
smooth plate channel. Figure 1 depicts the 2-D uniform staggered grid structure, in which the discrete pressure 
P is defined at cell centers, the discrete streamwise velocity u at centers of vertical cell walls, and the discrete 
normal velocity v at centers of horizontal cell walls. 
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u I V o p 


Figure 1. Staggered grid structure for 2-D channel flow. 


For planar channel flow, since the base flow is known (the Poiseuille flow), only the perturbation equations 
need to be solved. The governing equations described in section 2 can be simplified to two-dimensional equations. 
The computational domain is 


X € [0 ,mL x ], y e [0, 2], (67) 

where L x is the nondimensional wavelength of a 2-D T-S wave, and m is an integer. The base flow in this domain 
is the Poiseuille flow 



Mv) = i - (y- i) J . 

(68) 

and the resulting simplified 

governing equations can be written as follows: 


du 

<9(uo + u)u duv 1 d 2 u 5 2 u duo dP 

dx dy Re dx 2 dy 2 V dy + dx 

(69) 


dv d(tio 4- dvv 1 .d 2 v d 2 v. dP 

dt dx dy Re ' dx 2 dy 2 dy ' 

(70) 


du dv 

(71) 


dx + dy = °' 

Fourth-order central differences is applied here. Because of the use of staggered grids, 

we need the following 


two forms for the first derivatives of a generic function <^(x, y): 
• at x»: 


8<f) + + Hi+l,J - Hi-1, j + A _ ^ 

{ di k > ~ nZi + ° (Ax 

• at x i+i- 


(72) 


A 


+ 27^, +1|J - 27<ft t|J + 


+ <>( Ax*), 


24Ax 


(73) 
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The second derivative can be written as 


-&+ 2 , j + 

(ft?** ' 12A^ + ° (AX 

The formulas for the derivatives in the normal direction are similar. 


(74) 


The order of the finite difference approximation is reduced to second-order at the locations adjacent to the 
inflow/outflow boundaries to avoid very complicated one-sided high order finite difference approximations. Thus, 
the formulas are expressed as 


O'- = + 

• - ti±M. ~ 2< ^ 'id + tizhi ■ o(ax 2 \ 
~ Ax 3 +C»(AX). 


(75) 

(76) 


The normal derivatives were also changed to second-order for the points adjacent to the solid walls. The derivatives 
for the y-momentum equation (70) are approximated by 


■ **0^ + 0^,, 

(0k, = + AS). 


(77) 

(78) 


Since the no-slip boundary condition must be satisfied, we do not use the image points (located at y = — and 
y — 2 + ^ ), but derive the required second-order derivatives directly: 




2 _ M 
£ 2 




(79) 


( Qy2 ~ j > 


(80) 




v, 


J? 

7* 

ji « 

t*. 

aT 


L^r-j 








1 

J 

f ”>.iv 


(a) (b) 

Figure 2. Indices of the staggered grid. 




Before we express the resulting discrete system, it is necessary to clarify the special indices used for the 
staggered grid. The spacings are determined by 


Ax = 


mL x 


(ni - 2) ’ 


Ay = 7—7 


(nj — 2) ’ 


(81) 
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where ni and rxj are numbers of grid points in the streamwise and normal directions, respectively. Referring 
to the notation in Figure 2 (a), the resulting scheme may be confusing because of the existence of non-integer 
indices. Since three governing equations are discretized based on different reference points, it is possible to denote 
the discrete variables with three independent index sets (see Figure 2 (b)). Define 

= *i- 1. Vj ,=Vj\ ( 82 ) 

*». = *«. y,. = y J _« ; ( 83 ) 

this allows us to express the discrete system much more simply, with the computational domain now discretized 
using 

iu = (84) 

jv = 2 r",nj t (85) 

which means u 2 is located at inflow boundary and u\ j and Vi t j are ghost points located outside the computational 
domain. Similarly, is located on the lower solid wall and and V{ t i are ghost points. 


Since we have introduced the above three independent sets of indices, the discrete x-momentum equation is 
now based on (i Ul j), y-momentum equation on (i, j v ), and continuity equation on (i, j). Furthermore, denote 


UEE — u i* + 2,j* 
U NN ~ u t*,; + 2, 
UC ~ 

VEE — ^i+2,; v > 
VNN — v i t j* + 2> 
V C = 

PE = Pi+l.ji 
P N — Pi , j + 1 > 


U E — Ui^ + 1 J, 
UN — 

VE ~ 

VN = + 1 » 

PC — Pi,j> 
Ps — Pij - 1 ) 


U\v = 

Av = Pi-i.ii 
Pss — Pij-2, 


we obtain the following difference scheme written in general form: 


U\VW — 

U 55 = u i^,i-2i 


V 55 = 

Pww — Pi-2,jij 


A ^ e U £ £ F F Aty^w F F A^at^at// -I- Aatuat + 

A5U5 + Assess - Ac^c F APePe F APwPw H- APww Pww ~ APcPc = Sue, 
Bee^ee F Beve F F 5 ^%^ + ^nn^atn F -Bat^a^F 

i?sus F ^5s^ss — F BPnPn F BPsPs F BPssPss — BPcPc = Svc 1 
Cee u ee + Ce^e F CVujy — Cc^c F F D^^at F ^5^5 — Dc^c — •S'pc;* 


( 86 ) 

( 87 ) 

( 88 ) 

(89) 

(90) 

(91) 


Figures 3-5 depict the related variables for x-momentum, y-momentum, and continuity equations, respectively. 
In these figures, certain variables are obtained by fourth-order interpolation, for example, in Figure 3, 


(92) 


V C = 

= [ 9 (u 1? F Wij\+i F Vt-ij, F Vi-ij.+ i) “ (^i— 2,i„ — 1 F v»-2,j,+2 F + v t+i,j*+ 2 )]/ 32 , ( 93 ) 


with truncation error 0( Ax 4 -1- Ay 4 ). Figure 6 (a) shows the relationship between Vi^j and 


For points adjacent to solid walls, the accuracy of the interpolation in the y-direction is reduced to second 
order to avoid to use the ghost points which are impossible to be obtained accurately. The interpolation formula 
in this situation is (see Figure 6 (b)) 

+ Vi,y.+ i + - ( v >-2,j, + Ui-2,>.+i + v i+l,j. + t 'i+i,iv+l)]/ 32 > 

with truncation error 0(Ax 4 + Ay 2 ). 


(94) 


LiUfLiUfMcCormick , Final Report 



Figure 3. Neighbor points for the discrete x-momentum equation. 



Figure 4. Neighbor points for the discrete y — momentum equation. 



Figure 5. Neighbor points for the discrete continuity equation. 
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Vi-2 

,J. + 2 

Vi - 1 

J* F 1 

t 

1 

V i+ 1,3* + 2 




1 

Vi„j 

t 


Vi-2 

,i* - 1 



^t-1 


1 

ViJ. 

V i+ 

1 -. --F v j 

Vi-2 

,j.+i Vi - 1 

1 

1 

to 

,j«+i v »Fi,i»Fi 

t t 




1 

1 


1 

V iuJ 

1 

t 



i i 


V i-2,j* V i-l,j v v ij* V *Fl,j» 

(b) 

Figure 6. Interpolation for vatau point (a) interior point, (b) point adjacent to the solid walls. 


The resulting formulas of the coefficients for the interior grid points are listed below: 
• x-momentum equation 


Aee — 
Ae — 
A w ~ 
Aww = 
Ann — 
An — 
As = 
Ass = 
Ac ~ 
AP e = 
AP W = 
APww = 


12ReAx 2 
4 


12Ax 


(use F 


-(tifi F uqs), 


ZReAz 2 3Az ' 

4 2 / V 

-F — — \uw F u 0 *y)> 


3fleAz 2 

1 


3Ax 


\2ReAx 2 12Ax 
1 1 

12tfeAy 2 + 12 A ~y Vn 
4 2 

”n, 


(%iv F uowvv), 


3i?eAy 2 3Ay 
4 2 
3fleAy 2 + 3Ay l 


1 1 

~~ 12-ReAj/ 2 ~ 12Ay t '"’ 

3 5 / 1 1 V 

2A* + 2.Re ( Ax 2 + Ay 2 '’ 


24Ax ’ 
27 

24 Ax’ 


24 Ax ’ 
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A P c 


Su c 


27 

24 Ax’ 


-4 ug, + ** 


n— 1 
C 


2At 


- 2t/ c vcj 


where yc is the normal coordinate of uc* 


• y- momentum equation 


• continuity equation 


Bee — 
5* = 
= 

5wVV = 
Bnn — 
Bn — 
B s = 

#55 = 

5c = 
5iV = 
5P 5 = 

5^55 = 

BP C = 

Sue = 


1 1 / V 

+ 7X7 ( u ee + «Oee), 


125eAx 3 12Ax 
4 2 / X 

-(U e 4- TiOe), 


35eAx 3 3Ax' 

1 1 


125eAx 3 12Ax 

1 1 

+ 


{y>ww “I" ^Ouitu)j 


125eAy 3 12Ay 
4 2 

35eAy 3 3 Ay Ni 
4 2 

+ X7 — v s, 


VNNi 


3ReAy 2 3Ay 


1 1 
12i2eAy 3 “ 12A y VSS ' 

3 5 1 1 X 

2At + 25e Ax 3 + Ay 3 '’ 
1 

24 Ay’ 

27 

24Ay 5 

1 

~ 24Ay ’ 

27 

24Ay ’ 

~4^c H- ^c 
2A< 


Cee 


C c 


Dnn 


D c 

Spc 


-Cw — — 


1 


24 Ax 1 


C e = 


27 

24Ax’ 


-D s = 


24 A y* 


D n 


0. 


27 

24 Ay’ 


(95) 


(96) 


(97) 


Though Spc is zero at the finest grid level, it will be nonzero at coarser grid levels, so we still keep it in our 
formulas. 
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It is easy to describe the coefficients for the second-order accuracy scheme, but for those points adjacent to 
the solid wall, the coefficients are special. Only those for the i-momentum equation are irregular, so they are 
listed below: 

• lower wall: 


• upper wall: 


Here, 


Ann = Ass — As — 0 , 

4 1 

An = 


Ac = 


3ReAy 2 
3 


2A y n ’ 
1 / 5 

+ -57L . •> + 


2 At Re v 2 Ax 2 Ay' 


j) + 


2Ay 




Ann = Ass 
4 

= 


An — 0, 

+ 1 ' 


3 Re Ay 2 2A y ‘ 

, 3 1 5 

A ° ~ 2At + 2Ax 2 


+ 


Ay 2 ' 2 Ay 


v’ n = 0.5(i»ij.+i + i>i-x,>.+i)i 

v' t = 0-5Ki, +«<-!,>.)• 


(98) 


(99) 


( 100 ) 

( 101 ) 


4.2 Discretization on a Stretched Staggered Grid 

For flow over a flat plate, the computational domain is a y-direction truncated semi-infinite domain. Since the 
subdomain that we are most interested in is very close to a flat plate (scaled by boundary layer thickness), local 
high resolution near the flat plate is required, and thus a y-direction stretched grid which becomes increasingly 
dense near the solid wall is employed (see Figure 7). 



Figure 7. y-direction stretched grid. 


To meet the above requirement, an analytical mapping in the wall-normal direction is employed, i.e., 


y(v) = 


yrnoxcry 

ymaz{j1maz V) 


( 102 ) 
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or its inverse mapping 


»?(y) = 


VmaxV(^ + 1/mag) 
ymax{<7 + y) 


(103) 


where y max is the height of the computational domain in the physical coordinate y, r/ max is the height of the 
computational domain in the computational coordinate rj, and cr is a constant used to adjust concentration of the 
grid points. Thus, we can use the mapping 


_ ymax^V 

WmaxC “I" 2/max (f/max 


(104) 


to transfer the stretched grid in the physical (x,y) plane to the uniform grids in the computational (£,17) plane. 
With this transformation, the variables in the computational plane are now defined as 


U = yrju, V = v, (105) 

and the resulting governing equations, if we represent them by using only 17, V and P in the computational (£,??) 
plane, are: 


17-momentum equation: 


dU ld(UU + 2U 0 V) d UoV + UVq + UV. 1 ( d 2 U 1 d 2 f U. »,I7„ SP A 

* + ,-h( — i + *;< * ) - Ej* ' + * + + 8£ = °' < 106 ) 


V-momentum equation: 


dV d(U 0 V + UV 0 + UV) d(2V 0 V + VV) 1 . d 2 V 
yi) ~dt + 57 + 5 ^ + 


dt 


9r) 


Re'n 


dC 


1 d 2 v dv s dp n 

y v dr, 2+yvr,vv dr, )+ d^-°’ (1 ° ?) 


Continuity equation: 


dU dV 

7T - ~ 
d( dr) 


(108) 


The required derivatives in the governing equations can be obtained analytically: 


Vn 


Vyy 


’OmazymaT&{p “1" 2/max) 
[“Hmax^ “I" 2/max (*?m ax v)r 
2 1)max&{p *4“ 2/max) 

2/ma*(<7 +2/) 3 


(109) 

( 110 ) 


Based on the above derivation, all terms are considered in the uniform grids of the computational (£,77) plane 
only, and, thus, a procedure similar to that described in the last section can be used to derive the finite difference 
formulas. 


Figure 8 depicts the staggered grid structure used for the spatial discretisation in the computational (£,y) 
plane. Note that the discrete pressure P h is defined at cell centers, the discrete velocity U h is defined at centers 
of vertical cell walls, and the discrete velocity V h is defined at centers of horizontal cell walls. Moreover, the 
discrete continuity equation is associated with cell centers, the discrete U — momentum equation is associated with 
U h i and the discrete V — momentum equation is associated with V h . 
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Figure 8. Staggered grid for flat plate flow in computational (£,77) plane. 


The resulting discretized equations, if denoted symbolically, can be written as 

AeeUee + A&Ue + A\yUw + AwwUwW A AnN^NN -f A^UsA- 
AsU S + Assess — AqUc + Cww P\VW 4“ CyyPw 4" CnPn — CcPc 
BeeVee + BeV e + BwV\v 4- BwwVww 4- BnnVnn 4- BnVn + 
BsVs 4- Bss^ss ~~ BcVc + DssPss 4- DsPs 4- DnPn — DqPc 
FeeUee 4- FeUe 4- FwUw - FqUc + GnnVnn + GnVn 4- GsV$ - GcVc 


s U} (in) 

Sv, (H2) 

S Pc . (113) 


Since the process is similar to that in the last section, we only give the coefficients for the U-momentum 
equation (Figure 9). For the interior grid points, the formulas can be obtained as 


Aee — 
Ae = 
Aw — 


+ ; ( Uee + 2 Uoee), 


12i£eJA£ 2 12A ty Vy 

4 2 


3i?e$A^ SAfo,, 


(1 U B A 2 U 0E )> 


+ irrz — {Uw +2 Uq w), 


3 Re' 0 At 2 ' 3A^ 

Aww = ~ ( u ww +2U 0 ww), 


Awn — 
An — 
As — 


12Re' 0 AC 12A 

1 1 


+ 


URelArpy^y^ ' \2A t]y v 
4 2 


'(I'nn 4" Vbnn) 




ZRe^Arfy^y^ ^Ar)y m 
4 2 


+ 


ZRelArpy^y^ 3 A rjy Vp 


(V n + Von) + 

iV.+Vo.)- 


1222eJ Arjyjj' 
2 VvyVwy 


3 Re q Ayy^ 6 
ty^vv-y 
3 ReZAr)y Vp ’ 
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Here, superscripts n and n — 1 are used to indicate previous time values and the superscript n + 1 is understood 
for the current time, which is dropped for convenience. The lower case subscripts denote the approximate. Note 
that certain values above are defined by fourth-order interpolation (as in the last section). 
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Figure 9. Neighbor points for U — momentum equation. 


At the adjacent of upper and lower boundaries (say j u = 3 and nj — 2) except for the inflow and outflow 
boundary points, the finite difference scheme in the redirection is reduced to second order, while fourth order ac- 
curacy is maintained in the ^-direction to preserve the phase accuracy. The resulting coefficients for [/-momentum 
equation then are: 
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with and V[ obtained by fourth order interpolation (see Figure 10): 

K = + Vi-ij.+i) - {Vi+ij.+i + V-_ 2>i . +1 ))/16, 

v; = (9(V;, Jv + Vi-uJ - (Vi +1)i . + Vi-a,j.))/16, (116) 

and similarly for Vq and . 



Figure 10. Interpolation for and V[ at the adjacent of upper and lower boundaries. 


At the lower boundary, 
modified for (115): 


since a no-slip boundary condition is used, the associated derivatives need to be 
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At the upper boundary, a far field condition is assumed, which implies that all the disturbance should have 
disappeared, so we can also let U = V = 0 there, which yields 


An = 0 , 


A 5 


4 

2Re^Ary i y v _ j y r>p 


-(v: + v’) 


2Re' 0 Ar]y ne 


_ 1 (V/ + Vq,) 1 / 5 4 *7 yy T 

2At 2A n y,^ ^ Rx' 0 K 2A? ^ A^J ^ 2Re‘ Q Ar 1 ' 
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4.3 Outflow Boundary Treatment 

Outflow boundary conditions have been the focus of study by many researchers. Taking advantage of the 
staggered grid leads to a fairly effective approach, which we now describe. First, a buffer domain technique 
developed by Streett &; Macaraeg (1989) is applied to our problem. Thus, a buffer domain is appended to the end 
of the original outflow boundary to smear all possible reflections from the buffered outflow boundary (see Figure 
ii). 
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t = 0 


£ — L original 


£ — Ltotal 



Figure 11. Extended computational domain. 


The problem in general is that the conventional buffer domain is too long (usually four to eight T-S wave- 
lengths), which greatly increases computation. Our goal is to maintain the accuracy in the original computational 
domain, but eliminate all the possible reflection waves in a very short buffer domain. To realize the above goal, 
the governing equations in the buffer domain should be parabolized to have strictly outgoing waves. Thus, a first 
buffer function b(() is a introduced here and applied to the streamwise viscous terms: 


d 2 U 


HO 


d 2 U 

Up' 


d 2 V 

dp 


HO 


d 2 V 

W 


(119) 


fc(£ ) is a monotonously decreasing function from 1 to 0 so that the upstream effects of the streamwise viscous terms 
will gradually disappear in the buffer domain. The essential feature here is that all damping mode disturbances 
at the buffered outflow boundary become zero. To understand this, we rewrite the first equation of (119) as 


HO 


d 2 U 

~dp 


d 2 u 



( 120 ) 


Since &(£) — ► 0 at the buffered outflow boundary, and accordingly — ► oo, then, we can also consider that 


the buffered outflow boundary is compressed from { = oo by the function y/b(£). Now clearly, if the disturbances 
are stable (damping modes), they will vanish at £ = oo. To treat growing or unstable modes, we need a second 
buffer function 6^ e (0 that reduces the Reynolds number in the buffer domain gradually to less than the critical 
(or subcritical) Reynolds number and makes all the perturbation modes become damping: 


1 . fc*.(Q 

Re Re 


( 121 ) 


Thus, the new modified governing equations in the computational (£, ry) plane become: 
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dV d(U 0 V + UV 0 + UV) ( d(2V 0 V + VV) 
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< L 

iZeJ * ' Vn S £ 3 


1 


av\ sp 


+ y, 3 r 7 3 +y ' r> '" a J + 


dr) dr) 
dU dV 
d( + dr) 


0, 

0. 


The buffer functions are chosen as follows: 
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(123) 

(124) 


(125) 


It is clear that the first function decreases from one to zero very rapidly as one moves from the original outflow 
boundary to the buffered outflow boundary. The second function increasing from 1 to c + 1 is a quadratic 
function that is continuously differentiable at the original outflow boundary. Note that the total effect of these 
buffer functions is that, toward the buffered outflow boundary: 




the momentum equations become increasingly convection dominated in the (— direction, while the equations 
generally become parabolic; and 


• the momentum equations become highly diffusion dominated in the * 7 — direction. 

This treatment makes the outgoing waves propagate outword without reflection in the (—direction, and any 
oscillation in the 77 — direction will be effectively smeared in the buffer domain. 


Finally, we need to specify the buffered outflow boundary conditions under these modified governing equations. 
The parabolic character of the above equations require only two boundary conditions. As mentioned before, we 
have the disturbances tending to zero at the buffered outflow boundary (which is actually located at ( = oo), so 
one condition is 


P = 0 . 


(126) 


This is a very important condition since the elliptic character of pressure has not been modified in our new 
governing equations. Any improper boundary condition for P may cause trouble. For the second condition, we 
use the traditional extrapolation method for V , i.e., 


d 2 v 

w 


= 0 . 


(127) 


Though this condition may not be so accurate, accuracy of solutions in the buffer domain is not so important, and 
the main concern is that there must be no reflection wave traveling back to the original computational domain. 


Referring to Figure 12 , condition (126) is imposed directly on P by defining 

P c = 0 . (128) 

This is an implicit condition. With it, the discrete continuity equation associated with Pc is then used to define 
U at the buffered boundary: 

Vn ~ V c 

Ub=Uc -JL (129) 
Condition (127) is imposed by determining V at ghost points just outside the boundary: 


V E = 2 V c -V w . 


(130) 
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Figure 12. Buffered outflow boundary points. 

Note that the above treatment is only suitable for the perturbation equations. For the original equations, 
although the idea is similar, the resulting boundary conditions are treated differently. 


4.4 Distributive Relaxation and Line-Distributive Relaxation 

Developing efficient relaxation methods for the resulting discrete governing equations is a very important task 
for this research. Taking advantage of the staggered grid, we apply the distributive relaxation scheme originally 
developed by Brandt (cf. Brandt 1984) to our problems. 

Though our scheme originates from Brandt’s pioneering work, many modifications have been made for effi- 
ciency and convenience. Taking the transformed equations in the computational (£,tj) plane as the governing 
equations in this study, this special relaxation scheme performs in the following two stages. The first stage of 
this scheme involves point Gauss-Seidel relaxation on the central f7-point ( Uc ) for the associated discrete 17- 
momentum equation, followed by an analogous relaxation for Vc . The second stage involves a correction to all of 
the variables associated with the cell. Figures 13 and 14 depict the form of the corrections used in the second stage 
for interior and boundary points, respectively. The objective of this stage is to change Uc> V c , and P c accordingly 
so that the associated momentum equations are still satisfied, and the new residual of the associated continuity 
equation is zero (More precisely, the residuals of the associated momentum equations with frozen coefficients are 
unchanged, although u, v and p are changed to satisfy the continuity equation). 


By substituting the respective U E and Uc by Ue + z and Uc — e, Vn and Vc by + 6 and Vc - 6 , and Pc by 
p c _|_ ^ i n to both U and V momentum equations and eliminating the original equations, we obtain the following 
correction equations based on one grid cell: 


A E c+A c e-Cc 7 = 0 » 
B e $ + B c 6 — Del — 0) 

Fe{Ue + c) — Fc{Uc ~ + Gn{Vn + £) — Gc(Vc — 6) — Sm ■ 

To make the solving process easier, define 

^ — Cc{Bn + B C ) 

Dc {Ae + Ac) 


(131) 

(132) 


Then 

_ F e Ue ~ F c Uc + G n V n - G c Vc - S M 
F e -f Fc^ + Oc 
(A e -F Ac)w + (Be + Be) c 
7 = c^TDi 6 ' 

t — u)6. 


( 133 ) 
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After this, we modify 

Ue + <i 

Uc~t, 

V n + 6, 

V c -S, 

Pc +7- ( 134 ) 



Figure 13. Distribution of the corrections for distributive relaxation at an interior point. 


U E *— 
U c «- 

v N - 

V c «- 
Pc - 



Figure 14. Distribution of the corrections for distributive relaxation at boundary points: 
(a) inflow; (b) far field; (c) solid wall; (d) solid wall and inflow corner. 


Relaxation is based on the following steps through out all the grid cells: 

1. Freezing P and V, perform point Gauss-Seidel relaxation on (111) to obtain a new U in the whole compu- 
tational domain. 
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2 . Freezing P and U , perform point Gauss-Seidel relaxation on (112) to obtain a new V in the whole compu- 
tational domain. 

3. Modify the velocities U Ci U Ei V C) V N and pressure P c according to (132)-(134). 

The above procedure should be performed in some order on all grid cells. In our code, we sweep along horizontal 
lines from left to right (downstream), starting at the solid wall and proceeding to consecutive horizontal lines 
up to the far field. After each sweep, the coefficients for the discrete equations (11 1)— (1 13) are updated and the 
process is continued. 

The above point distributive relaxation works well when the ratio of the grid sizes, Ax/ Ay is not very large. 
Unfortunately, we can not always guarantee our grid ratio to stay in some specific range, and so, a more efficient 
relaxation scheme is desired. Physically, the most important issue is mass conservation, so, the new approach is to 
make the continuity equation converge as fast as possible. To realize this, a so-called line-distributive relaxation 
scheme has been developed in this study. 
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Figure 15. Distribution of corrections in the computational (£, 17 ) plane. 


The only difference between this new approach and the old point relaxation scheme is clearly shown in Figure 
15, in which we distribute the corrections not just in one single grid cell, but along a whole column of ceils. Based 
on this setting, the correction equations now become 


(A* e + A* c )ej — C J c jj — 0, 
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1 4- (B 3 N 4- B 3 c )6j — BpSj- 1 — ^c7j = 0» 
(rt + *cto + (°tf + Gt,)Sj-G* N 6 J+1 -G i c 6 i - l = s{ „ 

(135) 

where 
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•^0 

O 
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’fc? 
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H 

to 

(136) 

By defining 

= € jV^> 

(137) 

we have 

Cc (-®ir + ~~ ^N^i+ 1 ~ 

3 ~ D* c (4 + A’ c )6, 

(138) 

This formula is still too complicated and costly. Recall that the problem is time-dependent and At must be small 
enough to maintain accuracy. Thus, we have 


B 3 C >> B 3 Ni 

A 3 C >> A 3 E j 

(139) 

so that 

„ „ C ’c 4 

Ll) ft! : r- . 

7% Ac 

(140) 


This is only a constant. By applying (140) to (135), we obtain the following tridiagonal system: 

~G j c 6j - 1 + [(F 3 e 4* F 3 c )uj + (G 3 n + G 3 c )]6j - G* N 6j + i = 

j — 2, 3, * * • , My — 1. (141) 


Let 


a j = i^E + 4K + G 3 n + G * c , 

= —G 3 N) 
c j — —G’ C ' 

j — 2, 3, • • • , Yij — 1. 

The tridiagonal system now reads 
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Noting that 

u)j > 0 , F E > 0 , Fq > 0 , ^ ^ 


we obtain 


(142) 


( 143 ) 


|o 2 | > \b 2 \ > 0, |a,| > |6j I + \cj\, Ky-il > |c B ;-i| > 0, 
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and thus the above system can be solved by LU decomposition without pivoting very rapidly. With 6j computed, 
we can update U 3 and V 3 by 

U^U’s + S^, 

U 3 c^U’ c - 6j u>j, 

J = 2, 3, • • • , rij — 1, (144) 

and 


j - 3,4, • 1 n J - 1. 

Finally, the pressure corrections 7 j are determined by 


(145) 


72 = 


A * + A K 3 6 2 , 


C c 


= 7j-i + 


j = 3,4, - 1, 


and P is updated via 


P’c-P’c + 7>, 


3 — 2 , 3 , * • * , n j 1 . 


(146) 


(147) 


The programming process can be described as follows: 


• Freezing P and V , perform line (or point) Gauss-Seidel relaxation on (111) to obtain new U in the whole 
computational domain. 

• Freezing P and 17, perform line (or point) Gauss-Seidel relaxation on (112) to obtain new V in the whole 
computational domain. 


• For each i = 2, 3, * • • , ni— 1 in turn, and for all j = 2, 3, * ■ * , nj — 1 at once, change , Ui^+ij , V|j v , V^+i, 
and Pij using the corrections 6j , and 7^ in the form depicted in Figure 15 so that the corresponding 
momentum equations, (111) -(112), are unchanged and the corresponding continuity equations, (113), are 
satisfied. 


4.5 Multigrid Methods 

Multigrid is well known for its efficiency in attenuating all the errors with different wave numbers by relaxation 
on different grid levels. For our specific applications, several different techniques are developed to make the 
convergence process faster. 

4.5.1 Full- Coarsening 

For simplicity of discussion, we consider only the two-grid case as depicted in Figure 16. We use a full 
approximation scheme (FAS) to accommodate nonlinearities. A two-level FAS algorithm for an equation of the 
form 

L k W k - f k 

may be described loosely as follows: 

i) relax on L k W k = f k , 

ii) solve L 2k W 2k = L 2k ll k W k + 7j| fc (/' 1 - L k W k ), 

iii) replace W k <- W k + I k k {W 2k - ll k W h ). 


( 148 ) 
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The notation we have introduced includes the difference operators L h and L 2h , the restriction operators /£* (for 
the approximation) and I^ h (for the residual), and the interpolation operator I^ h . 



Figure 16. Two-level staggered grid (full coarsening). 


We first perform distributive relaxation for (11 1)— ( 1 13), then calculate the residuals for z, y momentum and 
continuity equations ( R on the fine grid. Relaxation on the coarse grid requires restriction of the 
approximations and residuals from the fine grid to the coarse grid to provide initial guesses and right-hand sides, 
respectively. For restriction of approximations, we use the following stencils: 

*?(») : [ I ] , 

w ■■ ii i ] ’ 

I™(P) : [ \ \ } (149) 

.4 4 * 

For restriction of residuals, we use the following full-weighting stencils, which come from the so-called area law 
developed by Liu (1989) (see Figure 17): 

il h (Ru) : ft f ! 1 

il h (Rv) ■ | | 

.8 8 J 

3 k (*m) = [til (150) 

.4 4 J * 

After all these restrictions are done, we perform coarse grid relaxation. 
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Figure 17. Restriction for 12*. 

The next step is to interpolate the corrections from the coarse grid to the fine grid: 

u h <- + /2 k (u)(u Jh - lf(u)« h ), 

+ v)v k ), 

ph pH + I$ h (p){p 2h - ll h (P)P k ). (151) 

Here we use bilinear interpolation, given by the following stencils: 

&( A«) : [ f 

L 4 J ’ 

4(Av) : [ | i ] > 

&(AP) : [ J J I (152) 

.16 16 J * 

Figure 18 depicts the relations between Au 2 * and Au*, Av 2h and Av*, and A P 2h and A P h corresponding to 
the above stencils. 



Figure 18. Bilinear interpolation for u, v and P . 


The FAS V-cycle based on the ingredients described here performed well in all of the cases we have considered 
for 2-D channel problems on this study. Typical convergence factors per V(2,2) cycle investigated were about 0.3 
for the fully-implicit scheme. 
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• p h ^ u h \ v h o P 2h — ► u 2h | v 2h 

Figure 19. Two-level staggered grid (normal direction semi-coarsening). 


4.5.2 Semi-Coarsening 

The full coarsening approach works well when the grid ratio Ax/ Ay is not too large. For a flat plate, because 
of the use of highly stretched grids, in some regions, this ratio may increase an order of magnitude or more. 
Under such circumstances, efficiency of the full- coarsening method described in the above subsection degenerates. 
To avoid this inefficiency, we use a semi-coarsening multigrid algorithm based on the modified point distributive 
relaxation described in the previous section. To accommodate the high degree of anisotropy caused by the dense 
concentration of grid points near the solid wall, we use a semi-coarsening scheme that involves grids that are 
coarsened in the vertical direction only. See Figure 19, which depicts the location of the velocity and pressure 
nodes for grid h and 2 h. Note the advantage here in treating the vertical boundaries by semi- coarsening since 
they coincide on all levels. 

The relaxation process is the same as described in the last subsection, but the stencils of restriction and 
interpolation are changed to follows: 


ll k (U) : [ 1 ] , 

■ 0 ' 

Il h (V): 1 , 

_ 0 _ 

I 2 h h (P) : [|], (153) 

This means that U 2h and P 2K depend only on their corresponding lower and upper nearest neighbors on the finer 
level grid, and V 2h is just V h evaluated at the same point (see Figure 19). For restriction of residuals, we use the 
following full-weighting stencils for semi-coarsening grids: 


il h (Ru ) : 

I 2 H h (Rv) : 



i 2 K h (RM ) : 


(154) 
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Here, R U} Rv and R M are the residuals for the i-momentum, y-momentum and continuity equations, respec- 
tively. 


We use bilinear interpolation for the variables, given by 

I 3 \(AU) : 


I 3 \(AV) : 


I 3 \(AP) : 
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3 

! 
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and 


and 
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!]■ 


(155) 


Let A be used to denote corrections, so that 

A P h = I k h {P 3k - ll hpk ) 

denotes the fine-grid correction to pressure, for example. Figure 20 illustrates these interpolation processes we 
use. For example, the corrections 


are determined as follows: 
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Figure 20. Bilinear interpolation for corrections to U, V and P. 


By using semi-coarsening multigrid, the convergence rate becomes much better the flat plate problem. An 
average convergence factor investigated under proper time step is about 0.2 for each V (2, 2) cycle. 
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5 Numerical Methods for Three-Dimensional Problem 


5.1 Discretization of the Navier-Stokes Equations on a Uniform Staggered Grid 

The perturbation form of the three-dimensional Navier-Stokes equations that govern the behavior of the 
fluctuating parts of the flow fields in a 3-D planar channel can be written as 


du Ouqu duu duv du$ duw 1 . d 2 u d 2 u ^ _ n 

- + — + — + ^ + 'jf + ^-telte + W + M )+ aZ- Q ' 


at 


dx 

dv 


dx 


dz 

dvw 


dvuo dvu 8vv 

dt + dx dx + dy dz 

dw dwu dwu 0 dwv dww 

dt ^ dx ^ dx dy dz 


_ j_ a** ^ = 0 

Re c?x 2 + dy 2 dz 2 dy 

1 d 2 w d 2 w d 2 w dP _ 

'te'dx* + dy i ' + 'd^’ + ~ ’ 

du dv dw 

+ tz + — = °- 




(156) 

(157) 

(158) 

(159) 


Here, u : v , w , and P denote the fluctuating parts of the flow field, and u 0 is the steady undisturbed Poiseuille 
flow. 



For temporal channel problems, periodic boundary conditions are assumed in both the x- and z- directions, 
while no-slip boundary conditions are imposed on the solid walls. Let L x and L t be the nondimensional wave- 
lengths of one T-S wave in the x- and 2 -direction, respectively. The nondimensional height of the channel is set 
to 2. Referring to Figure 21, the boundary conditions may then be written as 

U(&, 0,2,t) = V(XyOyZyt) = W(Xy0y Zyt) = 0, 

u(x, 2, Zyt) = v(x,2 yZyt) = w(Xy2yZyt) = 0, 

= «(* + L X yy } Zyt)y 
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v(x,y,z,t) = 
t o(x,y,z,t) = 
P{x,y,z,t) = 
u(x,y,z,t) = 
v(x,y,z,t) - 
w(x,y, z,t) = 
P(x,y,z,t) = 


v(x + L t ,y,z,t), 
w{x + L x ,y,z,t), 
P{x + L X) y,z,t), 
u(x,y,z + L z ,t), 
v(x,y,z + L z ,t), 
w(x,y, z + L z ,t), 
P(x,y,z + L z ,t). 


(160) 


A uniform staggered grid is used for the discretisation (Figure 22). Values of P are defined at its cell centers, 
u at the centers of the cell surfaces parallel to the (y, s)-plane, v at the centers of the cell surfaces parallel to the 
(x,z)-plane, and w at the centers of the cell surfaces parallel to the (z,y)-plane. 


In the same way as in Section 4, second-order backward Euler differences in time and fourth-order central 
differences in space are used to discretize system (156)- (159). 


y 


V 


u 


X 



Figure 22. Staggered grid structure. 


Using methods analogous to that of section 4, we can write the resulting scheme in the foil owing general form: 


Aeb^be + Aeue + Awu w 4- Aww^ww 4- Annu^n 4- Anun 4“ A5U5 4- Assttss-f 
Appupp -f- Apup 4- Abvb 4- Abb u bb — Ac^c 4- Dww Pww + DwPw 4- DePb ~ DcPc 
Bee^ee 4- Be^e 4- B\yvw 4- Bww^ww 4- Bnn^nn 4- BnVn 4- Bsv§ 4- Bss^ss 4- 
Bppvpp 4- Bpvp 4“ Bbvb 4- Bbb^bb — Bcvc 4- EssPss 4- EsPs 4- EsPn — EcPc 
Cee^ee 4- Cewe + Cwu>w 4- Cww^ww 4- Cn^wsn 4- 4- Csws + Css^ss 4- 

Cppwpp -}- Cpwp 4- Cbwb 4“ Cbb w bb — Cc^c 4- FbbBbb 4- EbPb 4- FpPp — FcPc 
DUeeuee 4“ DUe^e 4- DUwuw — DUquc + DVjtn v nn 4- DVnvn+ 
DV$vs — DVcvc 4“ DWppwpp 4- DWpwp 4- DWbwb ~ DWc w c 


S U9 (161) 
S V1 (162) 
(163) 
Sm,{ 164) 
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where the subscript C refers to the central point of the control volume, E to its nearest neighbor to the east, W 
to the west, N to the north, S to the south, B to the back, F to the front, EE is nearest point to E from the 
east, WW the nearest point to W from the west, and similarly for NN, 55, BB and FF. As an illustration of 
the notation we use, the relevant symbols for the discrete x-momentum equation are depicted in Figure 23. 


Though it is straightforward to list all of the coefficients and source terms in the above discrete system, we 
restrict ourselves to one set of the coefficients for equation (161) as an example. The coefficients and source term 
for the interior points of equation (161) with respect to uc now are listed below: 


Abe 

Ab 

Aw 

Aww 

Ann 

An 

As 

Ass 

App 

Af 

Ab 

Abb 

Ac 

Dp 

Dw 

Dww 

S u 


1 

” 12 -Re Ax 2 " 
4 

SReAx 2 “ 
4 

3fleAx 2 + 

1 

“ 12.ReAx 2 

1 

“ 12KeAy 2 
4 

3iieAy 2 ~ 


12 Ax 


(USE + VO Ee)i 


3Ax 

2 

3Ax 


[ve + uo £), 

(% + Vow)y 

i i 

(uww + Uqww), 


+ 


12Ax 

1 

12A y X 


3A y 


4 2 

3ReAy* + 3Ay 1 '*’ 

1 1 

_ 12iieAy J 12Ay 1 '"’ 

1 1 

~12Re&z* + 12Ay'" // ’ 

4 2 

3J7eAz 2 “ 3Al Wf ' 

4 2 

3ReAz* + 3A 

1 1 

“l^ReAx 2 12Az Wib ’ 

3 5 1 1 1 V 

2A t + 2JZe Ax 2 + Ay 2 + Ax 2 ^ 1 


24Ax’ 


D c 


27 

24Ax ’ 


1 

_ 24Ax’ 

— 4ug -j- Uq 

~~ 2A t 


2 ycv c - 


(165) 


Here, superscripts n and n - 1 are used to indicate values at previous time steps (superscript n + 1, indicating 
current time, is dropped for convenience), yc is the vertical coordinate of the node corresponding to uc- Lower 
case subscripts denote the approximate values of the v and w at points where the values of u with associated 
capital case subscript are located (see Figure 23), and should be obtained by at least 4th-order interpolation (see 
Section 4). 
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z x y x 

(•) 0 >) 


Figure 23. Neighbor points for x— momentum equation: (a) cross-section in (x, y)- plane, 
z-axis orthogonal to view, (b) cross-section in (e, z)-piane, y-axis orthogonal to view. 

For spatial simulation, the discretization procedure is the same, but the inflow and outflow boundary conditions 
should be treated very carefully to obtain reasonable results. This will be discussed in Section 5.4. 

5.2 Discretization of the Navier-Stokes Equations on a Nonuniform Staggered Grid 

A nonuniform grid must be employed for the problem with general geometry, or even for the boundary layer 
problems on a flat plate. Since transition is very sensitive to numerical error, the grid transformation, or more 
precisely, the Jacobian coefficients used in the governing equations, must be very accurate. It is believed that 
at least sixth-order accurate Jacobian coefficients should be used to maintain fourth-order accuracy where a 
numerical mapping is used. To avoid this complexity, we instead use the analytical mapping that maintains full 
accuracy. Since no orthogonality of the grid is required for the governing equations (see Section 3), we are able 
to use a simpler grid. 

5.2.1 Analytical Mapping 

First, we consider a very special but relatively simple mapping 

x = (, 

y = y({.*7 ,C), 

z = C, (166) 


or its inverse transformation 


t = z, 

v = v{x,y,z) 

( = z. 


( 167 ) 
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Under this special mapping, according to Section 3, we obtain 

J = T) y , 

t. = C* = 1 , 

6 =C* =0. (168) 

The remaining task now is to find a proper analytical mapping rj = rf(x ) y i z). Assume we are considering a 
flat-plate-type problem, so only one solid wall is located in the computational domain, the shape of which can be 
expressed as 


y 0 = f(x,z). 


(169) 


First, a rectangulation of the domain should be obtained. Letting y* be the associated intermediate variable, the 
rectangulation mapping can be written 


t _ (y f{ x i z ))ymax 

Vmax “* /(* i z ) 


(170) 


Under this map, the transferred domain becomes rectangular. Next, an analytically stretched mapping is employed 
to transfer the stretched grid in the intermediate (x,y',z) space to a uniform grid in the computational (£,77,C) 
space. In our study, this mapping is chosen as 


. _ fa + Vmax)y f 
^ cr 4* y 7 

The resulting mapping from the physical to the uniform computational grid is then just 

_ 3/mai( J 4 Umax )(l/ ~ f{ x t ^)) 
ymax{° + V) ~ f{*yZ){<J + y rnax y 

and its inverse mapping is 

_ wivmax - f(x> z)) 4~ ymaxf(v> z)(<T 4~ ymax ~ ??) 

1 lmax{& 4~ ymax V) 


(171) 


(172) 


(173) 


where y ma3 . is the maximum height of the computational domain in both the physical and computational space, 
and a is the parameter used for adjusting the density of grids near the lower solid wall. 


(173) expresses the complete map, which transfers the uniform grid in the (£,rj, £) space to the general grid in 
the (as, y, z) space that matches the solid wall. The useful derivatives that are required in the governing equations 
are now expressed below: 


Vx 

Vy 

Vx 

Vxx 

Vyy 

Vzx 


ymax J 4“ ymax)(l/ 2/max) 

\Vmax {p' + 2/) “ f{ x > ^X* 7 4" ymax)] 2 

ymax cr ( cr 4" ym ax Ky max 
[ymax H" y) — f{ x > ■ Z )( cr 4~ ymax)] 2 
ymaxfx®{p 4 ymax)(y ymax) 

[ymax{<7 + y) - f(x , Z)(a + y m ax)] 2 * 

l fxx [ymax( cr 4" y) _ f{ x i z)((T + ymax)] T 2 / 2 (^ + ymax ) 


— ymai ff (^ 4 2/max)(y 1/max ) ' 
~^ , ymax a { (T “I" ymax ) (ymax — 


[ymax^ + y) - f{ x >z){o 4" l/max)]^ 


[Vmaxip + y) - f(x, z)(cr + y maI )] 3 ' 

. f „ + v) - /(*, z)(<r + Vmo*)] + 2/* (a- + y mat ) 


= ymax<7{cr + ymax)(y - y max ) : 


[y m «x{o + y) - /(*, z)(<T + ymax )] 3 


where f z - f xx - §~4, /, - §£, and /„ _ §4. 


(174) 

(175) 

(176) 

(177) 

(178) 

(179) 
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5.2.2 Governing Equations under Special Mapping 

Under the above mapping in (171), the governing equations derived in Section 3 can be simplified. First, the 
variable transformation between the u, v, and w originally defined in the physical (x, y, z) space and the U, V 
and W originally defined in the computational ((, *7, C) space is now represented as 


u — Urj yi 
w = Wrj yi 

v = V-Ur) x - Wri M) (180) 

or 

U = u/rj yi 

W = w/rfy, 

V = (urj, + vr)y + wt),)/riy. (181) 

Most importantly, the transferred Laplacian operator in the computational space is relatively simple: 

+ hi + + + + 2r> ‘drjd( + {v * x +Vvv + (182) 

As we described in Section 3, in order to make the numerical procedure relatively simple, we consider both 
(u, v , tu) and ( U , V , W) as unknowns, and combine three momentum equations and one continuity equation to- 
gether with (180) or (181) to generate a closed system. The resulting momentum and continuity equations can 
now be written as follows: 

• Original equation (for base flow) 


du . dull duV duW s , d d 1 

d( + dr] + d( ^ + ^dl +ri *dJ P Re AlU ~ ’ 


3v .dvU dvV dvW. dP 1 A 

dt + Vv{ ~dr + ^r + ~dT ] + ^ ~ Re iv - °* 


dw ,dwU dwV dwW. . d d . n 1 

^ + ^-df + ^r + ^r )+{rt ^ + dc )p - k Aiw = 0> 


du dV dW „ 
~dt + lhi + lK ~ ’ 


(183) 

(184) 

(185) 

(186) 


where u, v, w, P, U, V, W are all variables. (183)-(186) together with (180) contains 7 equations for 7 unknowns. 
• Perturbation equation (for disturbance) 


du . d[u(U + U 0 ) + uoU] 3[u(l^ + Vo) + uo V] 

Qt + Vy( + Tr) 

d[u(W + Wo) + uoW) s( d d 1 

+ ) + (- + r, x -)P - — Aru 


= 0, 


(187) 


dv d[v{u + Uo) + vqU] d[v(V + Vo) + t» 0 V] 
M+Vvi ^ ^ 

i d[v(W + W 0 ) + voW], t _ dP 1 A .. 
+ ac } + R~e AlV 


= 0 , 


(188) 
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dw f d[w{U 4 Uo) 4 W 0 U] 8[w(V 4 Vo) 4 u>oV] 

+ M d( + Yr, 


. ^ + ^0) + ^], , , 

+ g- c ) + (V, 


d_ 

dr] 


+ r ( )p ~w^ w = 0 ’ 


dU dV dW 

l% + lhi + lK 


(189) 

(190) 


where u,u, w , P, Z7, V } W are all fluctuating parts of the corresponding variables, and uo, vo> Uo,V 0i Wo 
represent the base flow. 


5.2.3 Discretisation 

In the computational ((, 17, £) space, the grids are uniform. Suppose u, u, w and U t V j W are defined in 
terms of a staggered grid in the computational space (see Figure 24). Here, the values of P are associated with 
its cell centers, u and U with centers of the cell surfaces parallel to the (rj, () plane, v and V with centers of the 
cell surfaces parallel to the (£, £) plane, and w and W with centers of the cell surfaces parallel to the (£, r 7) plane. 
Note that though Figures 22 and 24 looks similar, they are defined in a different space. 



Figure 24. Staggered grid structure in the computational (£,*7, C) space. 

Second-order backward Euler differences are used in the time direction, and fourth-order central differences 
are used in space. Details of such an approach can be found in section 4. We can write the discretised governing 
equations symbolically as follows: 

Aeevee 4 Ae^e 4 Aw^w 4 Aww^ww 4 Asn^nn 4 Aprupf 4 ^ 5^5 -1- Assvss+ 

Appupp 4 Apup + Abvb + Abb^bb ~ Ac^c 4 Dww P\vw 4 DyyPw 4 DePe — DcPc = S U} (191) 
Bee^ee 4- Beve 4 Bwvw 4- Bwwvww 4 Bnnvnn 4 BpfVj y 4 - P 5 V 5 4- Bssvss + 

Bppvpp 4- Bpvp 4 - Bbvb 4- Bbb v bb — Bcvc 4- EssPss 4 EsPs 4 BnPn ~ EqPc — (192) 

Cee w ee 4- Cpwp 4 Cw^w 4 Cwwwww 4 C^n^nn 4 C^wpf 4 Csws 4 Css^ss 4 
Cppwpp 4 Cpwp 4 Cbwb 4 Cbb^bb ~ Ccwc 4 FbbPbb 4 FbPb 4 FpPp — FcPc — , (193) 
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DUeeUee + DUeUe + DUwUw - DUcUc + DVnnVnn + DV nVn+ 

DV S V S - DV C V C + DW ff W ff + DW F W F + - DW C W C = 5 m , (194) 

U C = Vy c U C i (195) 

= r) yc W c , (196) 

vc = V c -rh.Uc-rh m Wc (197) 


As an illustration of the notation we use, relevant symbols for the discrete (-momentum equation are depicted in 
Figure 25. Though both the original and perturbation forms are used, we only give the formulas for the latter 
here since formulas for the original equations can be obtained from the perturbation form simply by setting the 
values of base flow to zero. The coefficients and source term for the interior points of the discrete (-momentum 
equation (191) associated with u c are given as follows: 


Aee — 
A f — 
A w = 
Aww — 
A^n — 
An — 
A s = 
Ass — 
A FF = 
A F — 
Aq = 
A B b = 
Ac — 
De — 
D w = 
Dww = 
S tt = 


1 + Jh£-(u BS + 2uo£g ), 

12ReA( 2 12 A(' EE ri yBE h 

_i Vm. i 2u °^ 

ReA( 2 3A£ ' Ub + ~ ’’ 


). 


VyE 

4 , 2 rjyc (n , 2u ot y 

teA£ 2 + 3A£ ^ + *7 vW 

__ 1 ^(u ww + 2u °^ ), 

L2/ieA( 2 12A( Vyww 


otc 


+ 


2i2e Ar; 2 
iQc 2r? yC 


fyc /t r \ \r \ 

l v nn ~r K Onn / t „ n . i 

12Ar; 12iieA»7 


+ Sf< v - +M >‘)- 2w 


3ReAr)' 

“C 1»C m , v \ , TC 

12A^ ( " + °" )+ 12/ZcAt,’ 


'.Re At] 
1 


J , + M) 

Re AC 2 + \2AC n + T 7 V// 


. _ *M (W > + ^ 

3AC V 7 Vyf 

+ ^£ W + B* 

3AC fyt 


), 

), 


1 


(Wfck + 


VyC 

leA( 2 12AC 
5 1 «C7 

*■ 2Be A^ 2 + Arj 2 + AC 2 


>7yM 

1 


). 


). 


24 A£’ 

D c = 


27 

24A£ ’ 


24 A£’ 


-44 + u r‘ . -P„„ + 8P n -8P, + P i4 

2At + ^ C 12Arj + 

/ Inn + SuojsrV^ — 8uosVi + ^0S5^i# 

T]y C { 


12Ar? 
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-UoFFWff + SuopWf - 8uo B W b + Uo bbWm. 

+ 12AC 

1 / -Ufnn + 8Uf„ - 8lif, + U(„ -U(n» + 8u <n ~ 8u<, + U(,, ^ 

25 + mc aJ >■ 


(198) 


Here, superscripts n and n - 1 are used to indicate values at previous time steps, and superscript n+ 1, which 
indicates the current time step, is dropped for convenience. Lower case subscripts denote the approximate values 
of the v and w at points where the associated values of u with capital subscript are located (Figure 25). Other 
symbols used in the above formulas are as follows: 


a 

= vl + ril + vl, 

7 

= Vxz + Vyy + V* 


du 

u f 

r* 

ii 


du 


iW 

ii 


XI 


(199) 


All function values that are required at other than the canonical locations are obtained by fourth-order 
interpolations in the computational space. For example (see Figure 26), 

V c = (9(Vc + Vn + Vtfw + V\v) — ( Vsww + Vnnww + Vsb + Vnne))/^, (200) 

P t — (9(Ps + Psw) - (Pse + Psww))/1$' (201) 


Art 

Vnn 

Y-+- 

U NN 


U BB 


AC 

w tt , 

~*~U F p 



V n 

Pn 



w t 





h- ► 

“N 


up 

U\VW 

v c 

UW 0 _] 

_u c _ 

Ufl C 

U\VW 

W e 

U W 

_ 

U B 

o — 
Pww 

Pw 

— ► O — 

Pc 

Pb 


Pww 

Pw 

Pc 

Pb 


V. 

p. 

>-► 
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p.. 

USS 


C 


A( 


U BB 



! I 



0 *- ® 

C ( v C 

Figure 25. Neighbor points for (-momentum equation ( U are at the same points as u and are not shown here). 


The coefficients for the rj- and (- momentum equations are defined in an analogous way, and the discrete 
continuity equation is developed simply by applying central differences to each term. 

For the points adjacent to solid wall or the far field boundary, we change finite differences in the rj-direction to 
second order, and maintain fourth order in both the (- and directions. The method is the same as described 
in Section 4, and is thus omitted. 
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Figure 26. Neighbor points for fourth-order approximation for V c and P s . 



Figure 27. Extrapolation of ^ at grid point adjacent to solid wall. 


We depart from Section 4 in how we treat the gradient of pressure in the 17-direction at the point adjacent to 
the solid wall. Specifically, is estimated, using extrapolation. Letting i, Pi l , and P iu 1 be given 

by (201), P is fit by a second order polynomial role near the solid wall: 


P(r)) = art 2 + brj + c, 

where 


a 


( p i.,i-^,i + p i„0/(^n 2 ) I 

b 

= 


c 


(15P iw j-10P <%if +3P imit )/8 J 


This yields the extrapolated value 

,9P = 

' dr) '**•$ 2A»j 

A similar method is used to treat the point adjacent to the far-held boundary. 


( 202 ) 


(203) 


( 204 ) 
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5.3 Line-Distributive Relaxation 

The line-distributive relaxation scheme used here is basically the same as that described in Section 4. Though 
the specific governing equations are different, the general form is just (191)-(197). By using the predicted correc- 
tions distributed in Figure 28, for the Dirichlet boundary value problem we can describe the three-dimensional 
line-distributive relaxation as follows: 


• Freezing P, U , V, W, v, and w, perform point Gauss-Seidel relaxation on (191) over the entire computational 
domain to obtain a new u. 

• Freezing P, U , V , W, u, and w, perform point Gauss-Seidel relaxation on (192) over the entire computational 
domain to obtain a new v. 

• Freezing P, U , V , W, u, and v , perform point Gauss-Seidel relaxation on (193) over the entire computational 
domain to obtain a new w. 


9 Use the inverse transformation (195)-(197) to obtain new TJ , V, W . 

• For all j = 2 , 3 , • • • , nj - 1 at once: change U ijk) Ui + 1 j k ,V ijki W ijk , and W {j fc +i to satisfy the associated 
continuity equations, then update Pij k so that the new U , V, W and P as well as the associated transferred 
u, v , it; satisfy the three momentum equations. 


Since all of the u , v , w have been previously relaxed, and the U , V, W are updated, we may assume that 
equations (191)— (193) hold exactly. Let c, 6 , a and 7 represent the corrections for U, V , W and P, respectively. 
The correction equations for cube ijk corresponding to (191)— (193) are written as 


(4*C + ‘ ik + Ap^)e, - Dp 

Bp (6, - 6 i+ 1 ) - - 6 ■) - Epy, 

(Cp V y iih+t + Cpr,^-)a 3 - 

(DUp + DV'p)e } + ( DWp + DWpfa + DVg'fo - 6 j+1 ) - DVg'tfj-, - 6 ,) 

j =■ 2 , 3 , • • • , rij — 1 . 


0, (205) 

0 , (206) 

0, (207) 

SP, (208) 


This system has 4(n ; - 2) equations and 4 (n ; - 2) variables. Unfortunately, coupling between the correction 
variables makes the problem somewhat complicated. To develop a simpler approximate system, define 


with fixed i and k. Then 








V 


= 77 


Dp (Bp ± Bp)6, - BpS j+1 - BpS^ 
Bp (Apr,P'’ h + Apr,;^ 

Fp (Bp + Bp)S 3 - BpS j + 1 - BpS,., 


rij* 


(CpV? iik+ ' + CprX^b, 


(209) 

(210) 


From (198), we find that Ap ~ ~ for high Re and small At, which is much larger than Ap . Similarly, 

B P ~ jii and C P ~ jfa- These y ields 


u xi 




X1 E> c nP“ 

H 


E’^' 


At; 

A£t£* k t£*‘ ’ 

At; 

ACVy iih vV ih ' 


( 211 ) 


( 212 ) 
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Here, the superscript and to* denote the coordinate values of the corresponding discrete u 

points. 


Vjnjk = 0 


UiGk — ^6 

1 

Pi 6k + 76 

— O — 

V%6k + 6& — &6 

i 

ITii+1 6 ft + *6 

UiSk ~ *5 

P%Sh 
Vish 4 

+ 75 

— 

6<-6s 

tfii+1 5 * + c 5 

UiAk ~ ^4 

Piik 

+ 74 

> 

<$3 - 64 

^i +1 4 k + ^4 

U% 3 k — ^ 3 

Pi3 \ 

v i3k + 

+ 73 

) — 
62 — 63 

t^i+1 3 * + C 3 

Ui2k ~ € 2 

Pi2k 

-*■ C 

+ 73 

) — 
t 

Ui + I 2k + € 2 



Vi,* = 0 



Figure 28. Distribution of corrections in the (£ , rj) plane. 


With the above approximations, u) x j and ufjj can be treated as known parameters, so equation (208) < 
written in terms of only the unknowns 6j i 

\{DU'> k + DU'J k )cj Z} + (DVjf k + DV' c ik ) + (DW' F 3k + DW' c ik )w tj )6j - DV l c ik 6^ x - DVjf k 6 i+l = 

Let 

= (DU'J k + DU'J k )u xj + (DV$ k + DVy k ) + (DW' F jk + DW' c jk )w Mj , 
b 3 = -Dv' N > k , 
c, = -Dvy k , 

3 — 2| 3, * • • , n, — 1. 


. v, w 


be 


(213) 


(214) 

(215) 

(216) 
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Then we obtain the tridiagonal system 


CL2 63 

C3 <*3 


Cnj-2 Onj-2 ^nj-2 &nj- 2 Sm* 

C nj -i Onj-1 J . &nj- 1 J Y 1 * 

Thus, Sj , j = 2,3,*--, rij — 1 can be determined very efficiently. The other velocity corrections are given by 



— 

Uxjtj, 



j = 2, 3, • • • , rij — 1. 

on all cells 

in the i, k y— line as follows: 

U*+ 1 i k 



u ijk 

<- 

U ijk - 

\y'J *+i 

<_ 

W ij k+1 + a j, 

W i}k 


W' ]k - <Tj, 

3 = 2, 3, • • • , rij — 1, 

yij* 

<- 

i— t 

* i* 

1 c . 
t-H ♦ 

1 

+ *T 
II 


Finally, the pressure corrections jj are determined as follows: 


fa 4?* + B« k + 

^ 

%A¥ + Bii k + %C'S k 

(£+£+£)c ifc 

i = 3 , * • • , n j — 1 . 


P is then updated via 


Pijk <— P%jk + 7 > » ( 222 ) 

J — 2 , 3 , • • * , rij — 1. 

For large time step, the above method for pressure correction may cause some trouble because the assumption 

Ac >> Ae % Be >> Bffy Cc >> Cp 

is not valid. For such cases, we use the above method to obtain an initial guess for P, then solve an induced 
Poisson type problem for P with Neumann boundary conditions: 

d 2 P d , <9P, „ 

w + = ” 


Here, S p is the source term for the pressure equation. Only a few relaxations are needed to obtain a relatively 
good approximation to P, thus the increases cost is very limited. 
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5.4 Multigrid Methods 

Basically, the multigrid strategy for three-dimensional problems is the same as that for two-dimensional 
problems, although there are a few technical differences. To simplify the discussion, we limit ourselves to two 
grids. Note that the full approximation scheme (FAS) must again be used in order to accommodate nonlinearities 
of the Navier-Stokes equations. 

For a generic equation of the form 

L h W h = /\ 

a two level FAS process may be described loosely as follows: 

i) relax on L h W h = /\ 

ii) solve L 2h W 2K = L 2h ll h W k + ll h (f K - L k W h ), 

iii) replace W h <- W h + I$ h (W 2h - ll h W h ). 

Here we have introduced the difference operators L h (fourth order for fine grid) and L 2h (second order for coarse 
grid), the restriction operators (for the approximation) and I (for the residual), and the interpolation 
operator I k h . The coarse grid levels use the lower-order difference operator, which is less expensive and adequate 
for correcting smooth fine grid errors. Two coarsening strategies are used here: full coarsening, which is suitable 
for grid ratios Ax/ Ay near one and semi-coarsening, where the grid is coarsened only in directions of strong 
influence of the grid operators. 



u 3h | v 2K O w * h • p2K 

Figure 29. Projection of two-level staggered grid on the (x,y)-plane (z-axis is orthogonal to view, 
and the actual u,v } w and P points in the computational domain are located in different (x, y)-planes). 
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5.4.1 Full- Coarsening 


Figure 29 depicts the projection of a two-level full-coarsening staggered grid structure in the (z, y) plane. 
Under full-coarsening, the restriction and interpolation operators are defined as follows. For restriction l£ h of 
u, u, and ui, we use the stencil 


i r 1 1 

4 1 1 


(224) 


which means, for example, that u 2h takes one quarter 


l 2K 


A + U Ujl 


+ 1 *1 


*i = 2t - 2, 


from each u fc around it located in the same ( y , z)-plane: 

+ u £>i*,+i + u *ii+i*,+i)> 

j i = 2 j - 2, k\ — 2k — 2. 


For the pressure P, the stencil l\ h is 


1 

8 


1 1 
1 1 

1 1 ’ 
1 1 


(225) 


which means the restricted P 2h will come from all eight P* points contained in a cell of which P 2h is located at 
the center: 


p2/t 


g( P ‘ijikt + -Pii>i + 1 *i + + l ti+1 + 

P ii + 1 >i*i ^M + l >i + l ^»i + l >i*i + l "h ^m + 1 >i + l *i + l)» 

h = 2i - 2, ji - 2j - 2, ibi = 2fc - 2. 



momentum equations 



Figure 30. Full weighting restriction. 

For restriction of residuals, l£ h , we use the following full- weighting stencils (see Figure 30), which is based on 
the so-called volume law : 


1 

16 


1 

8 


1 1 
1 1 

2 2 
2 2 

1 1 
1 1 

1 1 
1 1 

1 1 
1 1 


for the discrete momentum equations, and 


for the discrete continuity equation. 


(226) 
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Taking the residual, r„, of the z-momentum equation as an example, we use 

r u ijk J^( r ** *i—i ji*i *** *i— i ii+i *i **« ii*i+i r * *i — i ii+i *«+i 
+ 2r u *i>i*i + 2r u i,i, + 1 k, + 2r i »iii*i + l + 2r i iiji + 1 *i + l 
*i + i ii*i r u *1+1 ji + i *i *1 + 1 ii*i+i *1+1 ji+i *1 + 1)’ 

where, ii = 2i — 2, ji = 2j — 2, ki = 2fc — 2. 

Trilinear interpolation is used for the transfers applied to the corrections 

A W k = W iK - ll h W k . 


(227) 


From Figure 29, taking u as an example, we first interpolate u on the (y, z)- plane, then on the (*, y)-plane. For 
the (y, z)-plane where both u h and u lh are located (Figure 31), the stencil is of the bilinear form 


' 9 
¥ 

3 ' 
¥ 

for Au£ w ; 

’ 3 
¥ 

9 ’ 
¥ 

. 16 

16 . 


. 16 

16 . 

' 3 
¥ 

H 9 * 

for Au*„; 

' 1 
¥ 

«| 9 »l 

- 16 

16 J 


. 16 

16 J 


for Au* e ; 
for AuJ e . 


(228) 


Once the Au l have been determined in the (y, z)-plane, all other points can be determined from these values 
simply by linear interpolation in the (z,y)-plane. For v and tn, the procedures are analogous. 
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X z 

Figure 31. Bilinear interpolation for Au in the (y, ;z)-plane. 

Interpolation, of the pressure corrections is quite complicated. The stencils can be expressed in block 
matrix form. Referring to Figure 32, the stencil for the specific A P h point is 


r 37 9 _ 

¥ 1 

64 64 


_ 9 _ JL 

¥ s 

64 64 J 


(229) 


where the first block of the matrix denotes the weights obtained from the rear plane, and the second from the 
front plane. 
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rear plane 


Figure 32. Trilinear interpolation for pressure points. 


5.4.2 Semi- Coarsening 

A full-coarsening strategy is generally ineffective for problems that favor special coordinate directions (e.g., 
anisotropic problems). To overcome this limitation, we consider now a special combination of semi-coarsening 
and line-distributive relaxation. The basic idea is, to use line-distributive relaxation in one direction (say the 
y-direction) and coarsening only in the other two directions (x- and 2 -directions). A two-level staggered grid 
projected in the (x,y)-plane and (x, z)-plane is given in Figure 33. 


Full weighting restriction is still used here for transferring the residual from fine to coarse grid. The stencils 
can be expressed as follows: 


n h (R») ■■ 
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( 230 ) 


These stencils can be explained geometrically as shown in Figures 34-35. 
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to 



to 


Figure 34. Full-weighting restriction 


for (a) i-momentum equation, (b) z-momentum equation. 



Figure 35. Full-weighting restriction for y-momentum and continuity equations. 


For the restriction of variables, bilinear interpolation is used. Its stencils are 
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For semi-coarsening, the coarse to fine transfer operators are based on linear interpolation: 
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(231) 


( 232 ) 


The meaning of the above stencils is also shown in Figures 36-38. 
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6 Results of Two Dimensional Simulation 

6.1 Two-Dimensional Spatial Channel 

We first perform two-dimensional spatial simulation for a planar channel. Let Re = 5000 and w = 0.33017. 
For these parameters, the Orr-Sommerfeld solution gives a = 1.1557 -fiO. 01 06. The perturbation amplitude is set 
to € — 1.2 x 10~ 3 . The channel length is set to five Tollmien-Schlichting (T-S) wavelengths, and one wavelength 
buffer domain is added to the rear of the channel. Thus, the total length of the computational domain is six T-S 
wavelengths. The grid used in this work is 192 x 64. Therefore, the grid for each T-S wavelength is 32 x 64. The 
grid is anisotropic with Asc = 0.1699 Ay = 0.03125. In order to maintain time accuracy, one T-S period is 
divided into 300 time steps. A semi-coarsening multigrid FAS scheme coupled with point distributive relaxation 
scheme is applied here. 

To examine efficiency, a comparison of the convergence histories of multigrid and single-grid relaxation at a 
fixed time step is depicted in Figure 39. The convergence factor investigated in this case is less than 0.1 per 
V(2,2) cycle, which is much faster than the depicted performance of single-grid relaxation. 

The streamwise and the normal components, u and v, respectively, of the perturbation velocity after 12 T-S 
periods are compared with the theoretical solution given by the linear theory at a vertical position close to the 
lower wall (y = 0.875 for v and y = 0.881 for u); see Figures 40(a) and 40(b), respectively. The excellent agreement 
in both amplitude and phase between our computational results and the theoretical solution was observed in the 
physical domain. The relative L 2 norms for both u and v (defined as and respectively, where u 

and v are numerical solutions and u and v are theoretical solutions) in the physical domain are 0.063 and 0.061, 
respectively, for the 192 x 64 grid. 

Figures 41 and 42 give the distribution of the respective disturbance vorticity and streamfunction at different 
times, showing that no order-one reflections were observed in the physical domain. The results are, of course, 
poor near the buffered outflow boundary, but they are located in the buffer domain and are therefore neglected. 


convergence history of mass (during one time step) 



0 


200 


400 


600 


Figure 39. Convergence history of semi-coarsening multigrid method. 
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disturbance velocity distribution (u) 



disturbance velocity distribution (v) 



Figure 40. Comparison of the numerical and theoretical solution at t=12T. numerical solution, 

^ 0 linear stability solution, computational domain: 5 T-S wavelengths physical + 1 wavelength buffer, 
grids: 192 x 64, Re = 5000, e = 0.001, scheme: fourth-order finite difference, flow is from left to right. 
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Figure 41. Disturbance vorticity contours for the first 3 T-S wavelengths of the 
physical domain, contour intervals=10“ 3 , flow direction is from left to right. 
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t=6T 



t=9T 



t=12T 


Figure 42. Disturbance streamfunction contours for the first 3 T-S wavelengths of the 
physical domain, contour intervals=5 x 10~ 5 , flow direction is from left to right. 


6.2 Two-Dimensional Spatial Flat Plate 

Since the linear stability theory (LST) is only valid for the parallel base flow, to test the accuracy of the 
scheme developed in Section 4, we first assume a ‘parallel Blasius profile’ as the base flow to compare our results 
with linear stability theory at the linear instability stage, then transfer to the more realistic non-parallel Blasius 
base flow. 

6.2.1 Parallel Base Flow 

Since linear stability theory (LST) only applies to the case when the amplitude of the imposed periodic 
disturbance is sufficiently small and the base flow is parallel, we first assume that the base steady flow is parallel 
to the solid wall, i.e., 

uo(*» v) = u o(*o, y), v 0 (x } y) = 0, 

where uo(xo,y) is the Blasius similarity solution at inflow. Since the finite difference schemes easily simulate 
decreasing (stable) modes, we restrict our tests to increasing (unstable) modes. 

Let itej = 900 and Fr — 86 (a; 0.0774). The Orr-Sommerfeld solution provides an eigenvalue 

a = a R + ia/ = 0.2229 - tO. 00451, 

and the associated eigenfunctions ^,^7, and <£/, which depend only on y. 

In LST, the disturbances are assumed to be travelling waves. In the 2-D case, this yields 

u = ee _0(rX (^cos(aiiX — ujRt) — (fffsin^aRZ — o>j*t)), 
v = ee _a/x [<f> v R cos(a rx — ujRt) — <f> v jsin(aRX — /*£)). 


( 233 ) 
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The inflow boundary velocities can then be obtained by setting x = xq for u and x — xq — for v (xo is the 
x-coordinate at inflow). We assume that e is small. Figure 43 shows the complex eigenfunctions for both u and 
v components (<££, <£/, <t> v & and <(>}) in the physical (x,y)- plane, which are normalized by setting max{^} = 1. 

A moderate size computational domain is selected. In this case, the plate length (measured from the inflow 
location x = x 0 ) is set to eleven ToUmien-Schlichting (T-S) wavelengths, while the buffer domain is an additional 
single wavelength, making the length of the total computational domain twelve T-S wavelengths. Here we choose 
y max = 75 and use a total computational grid with 362 x 50 grid points. The grid is highly anisotropic near the 
solid wall (Ax Ay). 


Since we use a fully-implicit scheme, the time step is restricted mainly by accuracy considerations and relax- 
ation efficiency. For our tests, we take 


2x 

320a; 


0.2537. 


The convergence rate of the semi-coarsening multigrid method coupled with point distributive relaxation 
scheme, which we generally found to be about 0.2 per V(2,2) cycle, is much better than the performance of 
single-grid relaxation, as shown in Figure 44. 


In order to compare our computational results with the linear theory solution, we first assume that the base 
flow is everywhere given by 

Uo((,r j ) = U o (0,T 1 ), V o ((,r,) = 0, and e = 10" 4 

(f/o(0,T7) is the Blasius similarity solution in the computational (£, tj) plane at £ = 0) and that the displacement 
thickness of the boundary layer is a constant, 6* o . Then Re*(= tfej), a, and w do not change along the streamwise 
direction. The streamwise and wall-normal velocity components, u and u, of the disturbance after 13 T-S periods 
(t=13T) are compared with the solutions obtained by the linear theory at a vertical position close to the solid 
wall ( y * = 1.3137 for u and y* = 1.2448 for u) in Figure 45. Excellent agreement in both amplitude and phase 
between the computational results of our fourth-order finite difference scheme and the solution obtained by LST 
is observed in the physical domain. 


To compare our numerical results with those obtained by LST more precisely, we check the profile of different 
Fourier modes by using the method mentioned in Section 6.3. Only the fundamental wave (k = 1) is obtained in 
this case. Figure 46 depicts the streamwise and wall-normal disturbance velocity profiles at x* = 439.2 and 608.4 
(x* = x/£* o ), and shows excellent agreement with those obtained by linear stability theory. On the other hand, 
the maximum streamwise and wall-normal amplitudes of the fundamental wave expressed by root-mean-square 
(rms) is also given in Figure 47. Very good agreement with LST is found again. 


To check for high order accuracy, two different grids (362 x 50 and 182 x 26) are used. For this kind of 
unsteady problem, it is quite difficult to assess grid convergence. Nevertheless, we attempt to do so by computing 
the relative L 2 error norm for both u and v given by 


11^112 = 






Il*ll» = 


/ £(«-«.)* 

5>« 2 ’ 


where u, v denote the numerical solution, and u e , v e denote the LST solution. The results given in Table 1 show 
that the accuracy of our scheme with this measure is about 0(h 3 S ). The maximum streamwise and wall-normal 
amplitudes of the fundamental wave for the 182 x 26 grid is also depicted in Figure 47 for comparison. 


Another important issue is the validity of the modified outflow boundary treatment for larger amplitude of 
disturbances. To test this, a series of disturbances with different amplitudes (10“ 4 to 0.1) are imposed. Note 
that, no order-one reflections (visual errors) were observed in the physical domain for all disturbances imposed 
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(see Figure 48). High frequency reflection maybe present, but do not affect the solution until sufficient energy 
has accrued in these modes. Note also that, for a large amplitude disturbance, linear stability theory is no longer 
valid and nonlinearity plays a significant role. 


grids 

ii-e.ii* 

l|E .|| 3 

182 x 26 

0.3315 

0.3243 

362 x 50 

0.0393 

0.0319 


Table 1 . Relative L 2 error norm for u and v after 13 T-S periods (calculated in physical domain). 

6.2.2 Non-Parallel Base Flow 

We next consider steady-state Blasius similarity solutions as the base flow. The local Re*, a, and u> now vary 
along the streamwise direction as the boundary layer thickness increases. We first use a 362 x 50 grid for the case 
of ReJ — 900 (ReJ = Uoq6* o /u is the Reynolds number at inflow) and Fr = 86 , which corresponds to a growing 
mode near the lower branch of the neutral curve. The increase of phase speed and amplitude along the streamwise 
direction is expected. In this context, our computational results are at least qualitatively correct. Figure 49 gives 
the computational results corresponding to Figure 45, and Figure 50 depicts the disturbance streamfunction and 
vorticity for this case at t = 13T. 

Overall validation of the 2-D code is made by carefully running a much larger case. To do this, an initially 
damping mode is imposed as the inflow boundary condition. The computational domain is selected to have a 
length of 32 T-S wavelengths in the streamwise direction (31 wavelengths physical domain and 1 wavelength 
buffer domain) and a height of — 75 (where = y/6^ and 6q is the displacement thickness of the boundary 
layer at inflow) in the wall-normal direction. The inflow data is given by: 

inflow Reynolds number ReJ = 688.315, 

frequency Fr = 86 (uj = 0.059195), 

amplitude of the disturbance e = 0.0025 rms. 

The linear theory provides eigenfunctions of disturbance for u and v, and the eigenvalue 

a = 0.174893 -f iO. 00501525. 

By using a 962 x 42 grid (30 grids for each T-S wavelength in the streamwise direction), we ran 33 T-S periods 
to obtain u(x,y, t) and u(x,y, t) at different streamwise locations. A Fourier transformation is again applied to 
compute the fundamental (T-S) wave ui, the mean-flow distortion uoi and the first harmonic wave 1 x 3 . 

Figures 51 depicts the maximum streamwise amplitudes of the fundamental wave Ui, the mean-flow distortion 
ii 0 , and the first harmonic 1 x 2 . Good agreement with DP-81 by Joslin et al. (cf. Joslin et al,. 1992; in DP-81, the 
grid size is 60/ x x 81, and l x refers to the number of T-S streamwise wavelengths included in the computational 
domain.) is observed. Figures 52 and 53 show the disturbance profiles of mean-flow distortion uo, u 0 and 
fundamental wave 1 x 1 , Vi with normal distance from the wall (defined as y* = y/S*) at the location Re* = 1519 
(Re* = UooS'/v refers to the local Reynolds number). Good agreement with DP-81 of Joslin et al. (1992) is 
again observed. 




LoglO(residual) 


60 


Liu } Liu, McCormick, Final Report 




Figure 43. Streamwise and wall-normal velocity eigenfunctions 
of an unstable mode ( Re J = 900) in flat plate flow. 
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Figure 44. Convergence history at a fixed time for multigrid and single-grid relaxation. 
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comparison of streamwise disturbance velocity 



Figure 45. Comparison of the numerical and theoretical velocity components near the solid wall (y* — 1.3137 
for u and y* = 1.2448 for u). Re* = 900, Fr = 86, parallel wall-bounded base flow assumption is used grids: 
362 x 50 (11 T-S wavelengths physical domain + 1 wavelength buffer domain). 
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10 -4 rms, x*=439.2 10~ 4 rms, x*=608.4 



10“ 4 rms, x # =439.2 1(T 4 rms, x*=6O0.4 


Figure 46. Comparison of the numerical and LST velocity profiles at x * = 439.2 and 608.4. 
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Maximum streamwise amplitude of fundamental wave u t 



X 


Maximum wall— normal amplitude of fundamental wave v t 



x* 


Figure 47. Maximum streamwise and wall-normal direction amplitudes of fundamental wave u\ and t»i with 
Re* 0 - 900, Fr = 86 and e = \/2 x 10~ 4 . 
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£=0.01 for u at j=0.8804878 e=0.01 for v at y=0.8193542 




€=0.03 for u at y=0.8804878 €=0.03 for v at y=0.8 193542 



e=0.05 for u at y=0.8804878 


€=0.05 for ▼ at y=0.8 193542 



Figure 48. Test of the efficiency of the modified outflow boundary treatment with different amplitude of distur- 
bances ( Re * = 900, t = 4T, grids: 66 x 50, 1 wavelength physical domain + 1 wavelength buffer). 
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comparison of wall-normal disturbance velocity f 



Figure 49. Amplitude of streamwise and wall-normal disturbance velocity components of 2-D disturbance near 
the solid wall ( y * = 1.3137 for u and y * = 1.2448 for u). Non-parallel Blasius base flow assumption is used, and 
the effect of non-parallel base flow is shown. 
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Figure 50. Instantaneous contour plots of (a) disturbance streamfunction and (b) disturbance vorticity at t ~ 13 T 
for iZeJ = 900, e = \/2 x 10 -4 , flow direction is from left to right. Non-parallel base flow is used, and only part 
of the computational domain is shown. 
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fundamental wave u, 




500 1000 1500 2000 



500 1000 1500 2000 

Re* 


Figure 51. Comparison of maximum streamwise amplitudes of fundamental wave ui, mean-flow distortion uo, 
and first harmonic uj with those obtained by DP-81 of Joslin et al. (1992) for Re 0 — 688.315, Fr 86 and 
e = 0.0025V2. 
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fundamental wave u t 



0.0 0.5 1.0 1.5 2.0 2.5 


mean— flow distortion u 0 



- 0.30 - 0.20 - 0.10 0.00 0.10 0.20 

10” 2 rms 

Figure 52. Comparison of streamwise disturbance profiles of fundamental wave Ui and mean-flow distortion 
uo against normal distance from wall between DP-81 of Joslin et al. (1992) and our approach (DNS+ MG) at 
Re* = 1519. 
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0.0 0.2 0.4 0.6 0.8 1.0 



0.000 0.002 0.004 

10~ 2 rms 

Figure 53. Comparison of wall-normal disturbance profiles of fundamental wave ui and mean-flow distortion 
uo against normal distance from wall between DP-81 of Joslin et al. (1992) and our approach (DNS+MG) at 
Re* = 1519. 
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7 Three-Dimensional Planar Channel Simulation 


7.1 Temporal Approach 

7.1.1 Linear Instability Stage 

Numerical simulations of transition in plane -parallel channel flow were performed for both the stable (Re = 
5000, a = 1.0, /3 = 1.005) and unstable (Re = 7500, a = 1.0, (3 = 0.3) cases. Figures 54 and 55 give the 
eigenfunctions for these two cases obtained from linear stability theory. 

Denote the fundamental wavelengths in the streamwise and spanwise directions by L x (— ^) and L t (= 
respectively. The imposed periodicity lengths are then L x and L M . The grid sise used was 34 x 130 x 18, which 
is anisotropic according to A z > Ax > Ay. 

To examine the efficiency of multigrid based on line versus point distributive relaxation, the convergence 
histories for the continuity equation at one time step are given in Figure 56, which shows that the line distributive 
method is much more efficient than the traditional point distributive method for this application, whether it is 
used alone or as a smoother. The multigrid method using line distributive relaxation has a convergence rate 
of less than 0.1 per V(2,2) cycle for appropriate At. Since the scheme is fully-implicit, it is much more stable 
than the explicit scheme. In all cases, only 125 time steps were needed for one T-S period, and on the NASA 
CRAY-YMP machine, only about 1300 seconds were needed for accurate computation of three T-S periods. 

The perturbation kinetic energy was used to check the accuracy of our scheme. The performance for 20 
T-S periods of the unstable mode and 3 T-S periods of the stable mode is illustrated in Figures 57 and 58. The 
respective increasing and decreasing rates of the perturbation energy show good agreement with the linear theory. 
Since we are more interested in the unstable mode, the streamwise, vertical, and spanwise components (u t v } and 
u>) of the perturbation velocity for the case Re = 7500 after 3 T-S periods are compared with the exact solution 
given by linear theory along selected x-~ lines. No significant phase errors are observed; see Figure 59. 

We also checked the relative La-norm error for the case Re = 7500 at selected x — y planes and different times 
during 20 T-S periods. The computation shows that during the first 1000 time steps (8 T-S periods), the errors 
change very slowly, but after that the errors become larger with increasing time. Since this is an unstable mode, 
the amplitude of the disturbance increases, so that the nonlinear effects may divert the numerical solution well 
away from that of the linear theory. The results are given in Table 2. 


step 

time 

IMa 

IIJMI* 


1 

0.59 

5.14e-2 

5.33e-2 

5.14e-2 

501 

98.51 

2.03e-2 

2.45e-2 

2.39e-2 

1001 

196.44 

3.33e-2 

3.49e-2 

3.66e-2 

1501 

294.37 

7.11e-2 

7.13e-2 

7.41e-2 

2001 

392.30 

1.10e-l 

l.lle-1 

1.13e-l 

2501 

490.23 

1.50e-l 

1.51e-l 

1.53e-l 


Table 2. Relative error Lj- norms for the z = 0.654 plane for u and v, and for the z — 0 plane for w. 
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Figure 54. Velocity eigenfunctions of a stable mode in plane- parallel channel flow (i£e = 5000, a = 1.0, 
(3 = 1.005). (a) streamwise velocity eigenfunctions, (b) real part of vertical velocity eigenfunction, 

(c) imaginary part of vertical velocity eigenfunction, (d) spanwise velocity eigenfunctions. 
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continuity equation, grid:34X130X18 



Figure 56. Comparison of convergence histories for relaxation only and for relaxation used with multigrid. 



Figure 57. Perturbation energy distribution in 3 T-S periods for Re = 5000. 



Figure 58. Perturbation energy distribution in 20 T-S periods for Re = 7500. 
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disturbance velocity distribution (u on y=l .5547. z=l 1.1265 line) 




0 Zn 

O o numeric*! folutioo theoretical solution 
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Figure 59. Comparison of the numerical and theoretical velocity at t — 3 T along 
(a) y = 1.5547, z = 0.6545 line for u and v, along y = 1.5547 , 2 = 0 line for u>, 
(b) y = 1.555,2 = 11.127 line for u and v, along y = 1.555,2 = 10.472 line for w. 
Re = 7500, a = 1.0, p = 0.3, e = 0.001, grid: 34 x 130 x 18. 
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7.1.2 Secondary Instability 

For secondary instability simulation, we use the so-called Benney-Lin-type (1960) initial conditions, 

V(x , y, *, 0) = Reoli^ivW + 4>{ K d \ (v)e iax+P ‘ + (vK“-" }, (234) 

as the initial guess for the velocity field. Here, the eigenfunctions fa d(y) <ta*±(y) correspond respectively 

to two-dimensional and three-dimensional eigensolution of the Orr-Sommerfeld equation at a specific Reynolds 
number, and the superscript k indicates the different velocity components, i.e., k = 1,2, and 3 denote the ti, v, 
and w components, respectively. 

The particular problem chosen here for study includes two cases: 

Re = 1500, a = 1.0, (3 = ±1.0; 

and 

Re = 7500, a = 1.0, (3 = ±0.5. 

Since the Re — 1500 case is linearly stable, a large amplitude of the disturbance is needed to obtain transition. 
Thus, we set the maximum value of the two-dimensional disturbance to 0.1H7o> and that of the three-dimensional 
disturbance pair to O.OSC^Oj where Uq is the centerline velocity of the undisturbed mean flow. For the Re — 7500 
case, which is linearly unstable, we set the maximum value of the two-dimensional disturbance to 0.03C/o, and 
that of the three-dimensional disturbance pair to O.OlI/o- 

We restrict the computational domain to 

*e [<>,—], ye [o,2], 

Ot 

Thus, for Re — 1500, the computational domain is 

0 < x < 2 t, 0 < y < 2, x < 2 < 3*. 

Total vorticity, which is defined by (cf., Biringen, 1987) 

n = ^/uf+u> J, 

is used as a metric to describe the behavior of the flow, where 

dw dv 
dy dz 

is the streamwise vorticity and 

_ dv du 
dx dy 

is the spanwise vorticity. Figure 60 gives the contours of the total disturbance vorticity at t = 15.8 in some 
selected (x, z)-planes. The process of A wave generation, is shown clearly in these figures. Figure 61 gives the 
3-D contours of the vorticity magnitude, showing clearly the A wave structure. 

The streamwise vorticity itself is another metric often used to analyze computational results. The streamwise 
vorticity contours in some selected (z, y)-planes at t = 15.8 and t = 29.1 are given in Figures 62 and 63. Although 
the grids and time are different, the results shown in Figure 62 are qualitatively comparable to those obtained by 
Zang et al. (1987). 

We also compared our results with the experiments (Nishioka et al. 1981) by using the equishear lines 
(which correspond to the approximate spanwise vorticity u M ) in the (x,y)-plane at the position of maximum 


T 3t. 

ze[ W\' W 
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Results from the computation that correspond to various stages of the laboratory flow between the lower wall 
(y = 0) and the channel center (y = 1) are presented in Figures 64-68. Figures 65-67 correspond to Figures 
4-6 of Nishioka et al.(1981). The Figure 64 contour plots correspond to the beginning of the computation, 
and the Figure 65 contour plots correspond to the one-spike stage (t=0.15 and t=7.99). Figure 66 shows the 
equishear lines at t = 15.8, corresponding to the three-spike stage. Due to the secondary instability manifested 
in the previous stage, breakdown of flow into smaller scales begins to appear in the computation. Contours of 
equishear at t — 23.6 are given in Figure 67, which correspond to the flve-spike stage of the experiment. Both 
results of numerical simulation and experiment show the intense shear layer developed in the wall region. The 
most significant feature of Figure 67 is the existence of distinct vortex structures in the wall region in both the 
laboratory flow and computation. Figure 68 gives the results obtained at t = 27.5 and t = 29.1, showing clearly 
that the flow scale has become much smaller, with large scale vortices broken down into many small vortices. 

Unlike the Re = 1500 case, the flow is linearly unstable for Re = 7500, so, even very small disturbance can 
induce transition. Similar results for this case are shown in Figures 69-71. Here, the computational domain is 
set to 


0 < x < 2?r, 0 < y < 2, < z < 6 t. 

Since the amplitude of the disturbance is much smaller than that for the case Re = 1500, the high-shear layer 
seems to spread more slowly at the beginning of the simulation. Laboratory experimental stages are not so 
clear in our selected time steps, but the qualitative behavior of the flow is still the same. The high-shear layer 
develops and then rolls up, finally resulting in the laminar flow breakdown. Because of the resolution limitation, 
we still cannot obtain the turbulent spot. Figure 69 depicts the contours of the total disturbance vorticity in the 
(x, z)-plane at t=47.3, showing A wave generation. Figure 70 depicts the contours of the streamwise disturbance 
vorticity, showing clearly the roll-up procedure. Finally, Figure 71 gives the equishear contour plots at t — 47.3. 
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Figure 60. Total disturbance vorticity contours at 0.1 intervals for t = 15.8, Re = 1500, a ~ 1, and f3 = ±1 on a 
34 x 66 x 50 grid. 
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Figure 61. 3-D contour plots of the vorticity structure (relative helicity) for the three-dimensional temporal 
channel a 34 x 66 x 50 grids. Re = 1500, a = 1, /? = ±1. 
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Figure 62. Streamwise disturbance vorticity contours at 0.20 intervals for t = 15.8, Re — 1500, a = 
/3 — ± 1 on a 34 x 66 x 50 grids; dashed lines indicate negative contours. 
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Figure 63. Streamwise disturbance vorticity contours at 0.50 intervals for t = 29.1, Re = 1500, a = 1, and 
= il on a 34 x 66 x 50 grid; dashed lines indicate negative contours. 



Figure 64. Contour plots of du/dy at t = 0.15 in the (x, y)-plane. Re = 1500, a = 1, and (3 = ±1. The contour 
interval is 0.2, flow direction is from left to right. 
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Figure 65. Contour plots of du/dy\ (a) one-spike stage, experimental results by Nishioka et al.(l98l); 

(b) computations at t = 7.99 in the (x, y)-plane, contour interval is 0.2. Flow direction is from left to right. 



Figure 66. Contour plots of du/dy: (a) three-spike stage, experimental results by Nishioka et al.(l98l); 
(b) computations at t = 15.8 in the (x,y)-plane, contour interval is 0.2. Flow direction is from left to right. 
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Figure 67. Contour plots of du/dy : (a) five-spike stage, experimental results by Nishioka et al.(198l); 
(b) computations at t = 23.6 in the (x,y)-plane, contour interval is 0.2. Flow direction is from left to right. 


Figure 68. Contour plots of du/dy at: (a) t = 27.56; (b) t = 29.1 on the central (z, y)-plane 
Flow direction is from left to right. Results show that the high-shear layer has broken 
and spread to the whole domain, and the scale of the flow has become smaller. 
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Figure 71. Contour plots of du/dy at t - 47.3. Re = 7500, a = 1, (3 = ±0.5. Flow direction is from left to right. 


7.2 Spatial Approach 

To accommodate the treatment of the inflow and outflow boundaries, we use a semi* coarsening multigrid 
approach where coarsening occurs only in the (y, 2 )-plane, which has the additional benefit that restriction and 
interpolation processes are simpler than those for full coarsening. 

7.2.1 Linear Instability Stage 

To test the accuracy of the scheme, we first restrict our problem to the linear instability stage. Letting 
Re = 7500, u) — 0.25, and /? = 0.3, then the least stable mode can be obtained as the eigenfunction of the 
Orr-Sommerfeld equation. We use this eigensolution as the inflow boundary condition. The amplitude of the 
disturbance is set to 

moz{^i}= 10" 4 
Vo 

By using a 182 x 66 x 18 grid, which includes a five T-S wavelengths physical domain and a one wavelength buffer 
domain, and dividing one T-S period into 200 time steps, we obtain exceptionally accurate results in the linear 
stage (see Figure 72). Note that the buffer has successfully damped all possible visible reflections from the outlet 
of the channel. To be more precise, we used the relative L 2 norm of the error 


e 2 (4>) = 



m 


where <f> denotes the approximate solution and <f> t denotes the exact solution. Results are given in Table 3. 
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norm 

E 2 (u) 

E 2 (v) 

E a (id) 

error 

0.033 

0.032 

0.041 


Table 3. Relative L 2 norms of the numerical solutions at t = IT 
(Re = 7500, and a 30 x 66 x 18 grid is used for each wavelength). 


7.2.2 Secondary Instability 

To study the secondary instability problem, we use the same problem formulation for which the inflow bound- 
ary condition is still of B-type, which includes a 2-D wave with amplitude 0.03 and an oblique wave pair (/? — ±0.3) 
with amplitude 0.01: 

V(0, y, z, t) = Real{<j>[ K J (y)e iut + (y) c -‘+*‘ + (y)e*“‘-^}, 

where <f>2d(y) and <f>3d± (y) correspond respectively to 2-D and 3-D eigensolutions of the Orr-Sommerfeld equation 
at Re = 7500, u = 0.25, /3 = 0 and ±0.3, respectively. 

Figure 73 gives three depictions of the equishear contours at different (x,y)-planes at t = 91.23, clearly showing 
the process of high-shear layer roll up. Figure 74 shows the total disturbance vorticity contours. 


distribution of u along y=0. 570 l,z=30. 7614 line 



Figure 72. Comparison of the numerical (<>) and theoretical ( — ) solutions at selected i— lines at t = 7T on a 
182 x 66 x 18 grid (5 T-S wavelengths physical domain ± 1 wavelength buffer), with Re — 7500, u; = 0.25, and 
P — 0.3. Flow direction is from left to right. 
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x= 16.3259 



Figure 73. Contour plots of du/dy at 0.30 interval for selected (x f y)-planes on a 162 x 66 x 50 grid at t — 91.23 
(7 T-S wavelengths physical domain -1- 1 wavelength buffer domain), with Re — 7500, u> — 0.25, (3 — ±0.3. Flow 
direction is from left to right. 


y=0.l25 




Figure 74. Total disturbance vorticity contours at 0.05 intervals for selected (x, z)-planes on a 162 x 66 x 50 grid 
(7 T-S wavelengths physical domain + 1 T-S wavelength buffer domain), with Re — 7500, w — 0.25, f3 — ±0.3, 
and t = 91.23. Flow direction is from left to right. 
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8 Three-Dimensional Flat Plate Simulation 


8.1 Linear Instability Stage 

We first perform the simulation for the small disturbance case to check the accuracy of the code. The 
computational domain is restricted to 


X £ [x 0 , Xq “1“ lx Ao] i V ^ 1 /mfli]) 


z € 




w m 


(235) 


where l x is the number of T-S wavelengths in the computational domain, and Ao is the T-S wavelength at inflow 
(the T-S wavelength A varies when the base flow is non-parallel). 


8.1.1 Parallel Base Flow 

In order to compare our results with those obtained from linear stability theory, a wall-bounded parallel base 
flow is again assumed, and a truncated computational domain with y max — 50 is used. Let ReJ — 900, Fr — 86, 
and the spanwise wavenumber (3 — 0.1. The LST gives the eigenvalue 

ot 3d = 0.2169 - iO. 00419. 

Computations are performed on a 170 x 42 x 18 grid (including a five T-S wavelengths physical domain and a 
one wavelength buffer domain), and each time period is divided into 250 time steps. By setting the amplitude 
of the disturbance e = 10~ 3 , and the stretch parameter a— 4.25 (note in (173), f(x,z) = 0 for the flat plate), a 
computational result is obtained which has good agreement with LST (see Figure 75). It can be founded that 
the DNS solution still has some difference when the flow propagates downstream, this is mainly caused by the 
too coarse grid (especially in the spanwise direction) we used. 


8.1.2 Non-Parallel Base Flow 

Second, a non-parallel Blasius base flow is applied, and the same grid and same computational domain are 
used. As expected, an increased growth rate is observed (see Figure 76). This is at least qualitatively comparable 
to linear theory, which exhibits a higher amplification factor with higher Reynolds number under the unstable 
region. 


8.2 Secondary Instability 

Finally, to demonstrate the real ability for simulating the secondary instability, a Benney-Lin type disturbance is 
imposed on the inflow: 

V(0,y, z,t) - Real{c 2 d4>^d {y) e%Ut + e 3d + ^3d+ (y)e twt+, ^ a + *3d_<^3d^ (v) e%ut 

where <f>2d{y) and < p 3 d ± {y) correspond respectively to 2-D and 3-D eigensolutions of the Orr-Sommerfeld equation 
with Re* 0 = 900, Fr = 86, (3 = 0 and ±0.1, respectively. The height of the computational domain is set to 
2/max = 75. A 202 x 50 x 34 grid (including a seven wavelengths physical domain and a one wavelength buffer) is 
used. The time step is set to 1/250 of the 3-D T-S wave period, and the amplitude is set to c 2 d = 0.03 for the 
2-D wave and 63^ = 0.01 for the 3-D wave. 

The code has been proved to be very efficient. About 30 words of memory is required per physical point, and 
about 40 ps of CPU is required on the CRAY-YMP per physical grid point per time step. According to estimates 
by Herbert (1991), direct simulation for the existing methods requires on the order of 100 megawords of memory 
and 1000 CRAY-YMP CPU hours for limited resolution. However, our approach only requires a total of about 10 
megawords of memory and 5 CRAY-YMP CPU hours to obtain approximately the three-spike stage of transition 
(Figure 77). 


86 


Liu, Liu, McCormick, Final Report 


Figure 77 gives two depictions of the spanwise vorticity contours on the central (z,y)— plane at i = 244.17 
and 270.97, which clearly show the process of high-shear layer roll up and the occurrence of spikes. Figure 78, 
showing part of the computational domain, gives the equishear ( —• ) contour plots at t = 270.97, showing good 
agreement with the experimental result obtained by Kovassnay et al. (1962). Figure 79 depicts the total vorticity 
contours at t = 244.17 and 270.97, which is also qualitatively comparable to the laboratory experiments obtained 
by Saric k Reed (1987) and Knapp k Roache(1968). 



303.9 336.7 373.5 406.2 443.0 477.8 




Figure 75. Comparison of the DNS and LST velocity components near the solid wall (y5 = 1.719, Zq = 32.84 for 
u; y5 = 1.624, *5 — 32.84 for u; y5 = 1.719, Zq — 31.41 for tu) on a 170 x 42 x 18 grid. Re 5 = 900, Fr = 86, 
€ — 10~ 3 , parallel wall-bounded base flow assumption is used. 
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physical domain buffer domain 


Figure 76. Streamwise, wall-normal, and spanwise disturbance velocity components of a 3-D disturbance near 
the solid wall (y£ = 1.719, = 32.84 for u; y* = 1.624, = 32.84 for v; = 1.719, z* = 31.41 for w) on a 170 

x 42 x 18 grid. .Reg = 900, Fr=86, e = 1 0 — 3 , non-parallel base flow is used, and the effect of non-parallel base 

flow is shown. 
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t=244.17,contour interval=0.025 



(a) 



(b) 


Figure 77. Spanwise vorticity contours on the central (a;, y )— plane of the computational domain at (a) t— 244.17, 
and (b) t=270.97. 



W 



Figure 78. High-shear layer at one spike stage of transition in boundary layer flow, (a) Numerical result in a 
selected part of the central (i, y)-plane, (b) experimental results obtained by Kovasznay et al. (1962). 
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Figure 79. Total vorticity contours on the selected (x,z )— plane at (a) t = 244.17, y£ = 0, (b) t — 270.97, yj = 0, 
Re$ = 900, grids: 202 x 50 x 34, contour interval is 0.035. (c) Aligned pattern of A-vortices (experimental results 
by Saric & Reed (1987)). Flow direction is from left to right. 
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9 Effect of Roughness Elements on Transition 

In this section, we investigate flow instability in the presence of one or multiple surface roughness elements in 
both 2-D and 3-D channels and flat plate boundary layers. Some of the computational results are compared with 
the experimental results, including that for a 2-D isolated large-scale roughness element, which is quantitatively 
compared with experimental results obtained by Dovgal & Koilov (1990), and that for 3-D single and multiple 
roughness elements, which are qualitatively compared with experiment given by Mochizuki (1961a,b) and Tadjfar 
et al (1993), respectively. For other cases, we just describe the results and give some qualitative explanation 
because of the lack of existing data in this new regime. 

9.1 Two-Dimensional Roughness Element for Channel Flow 

We first test the effect of isolated and multiple roughness elements for 2-D channel flow. The grid generation 
process was described in Section 5.2. For the 2-D case, the shape of the roughness element can be expressed as 
y 0 — /(x), which results a simpler mapping according to (172) and (173) and a nonorthogonal gird in the physical 
(x,y) plane. 

9.1.1 Isolated Roughness 

Let Re = 5000, which corresponds to a decreasing mode according to linear stability theory. To use our 3-D 
code for a 2-D channel case, we specify a 170 x 50 x 4 grid, including a five T-S wavelengths physical domain and 
a one wavelength buffer domain. The height of the channel is set to 2. The roughness surface is given by 

f(x) = 0.1h/cosh 2 (y/2(z - x 0 )), 

where xo is set to the 30th grid point in the streamwise direction. The stretch parameter o is set to 4 in the 
y-direction to induce more grid points in the roughness region. The grids are shown in Figure 80 which are slightly 
stretched and nonorthogonal to the roughness surface. Note that we use an analytic function f(x) as the surface 
function to maintain high accuracy in the coordinate transformation. 

Figure 81 displays the contour plots of the streamfunctions of the flow field without the imposed disturbances. 
A separation region appears after the roughness element. Using this flow field as the base flow, and setting the 
amplitude of the disturbance e — 0.0025 y/2, we find the disturbance becomes unstable in certain regions. Figure 
82 gives the instantaneous contour plots of the perturbation streamfunctions, showing increase of disturbance 
occurring in certain regions. Because of the lack of further excitation from the recirculation bubble, this increase 
cannot be maintained very far. A detailed analysis by Fourier transformation shows that the mean flow distortion, 
fundamental wave, and the first harmonic wave have abrupt increases past the roughness elements due to the 
effect of the recirculating bubble, which is in contrast to a smooth planar channel, that exhibits only a decreasing 
fundamental wave. The results are depicted in Figure 83. 

9.1.2 Multiple Roughness Elements 

To test the effect of the height of the roughness, as well as the effect of multiple elements, a 402 x 66 x 4 grid 
(including a nine wavelengths physical domain and a one wavelength buffer domain) is used. Here we set the 
height of roughness elements = 0.12, the roughness surface is described as 

m 

/(*) = ^0.12/coa>i 2 (2(z - *,)), 

1=1 

and o = 4.5. The first roughness element is at i = 82, which is two T-S wavelengths down from the inflow 
boundary. We first test the case with two roughness elements in this grid setting, with the second roughness 
located at i = 122. The results are shown in Figures 84-86. Since the region downstream is smooth, the 
disturbance becomes damping again after the flow is well past the roughness region, and the higher frequency 
waves damp even faster. 
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Finally, we placed seven roughness elements in the computational domain, starting from i — 82 and spaced 40 
grid points apart. The results, given in Figures 87-89, show very clearly that the disturbance is amplified after 
it passes each roughness element, apparently due to its receiving energy from the mean flow. 




Figure 81. Contour plots of the streamfunctions of the base flow with one 
roughness element. Re = 5000, k = 0.15.Flow direction is from left to right. 



Figure 82. Instantaneous contour plots of the perturbation streamfunctions for the one roughness 
element case. Re = 5000, k - 0.15, c = 0.0025\/2- Flow direction is from left to right. 
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Figure 83. Maximum amplitudes of fundamental wave ux, vi, mean-flow distortion uq, v 0} and first harmonic 
112 , V 2 for Re = 5000, k * = 0.15, and c = 0.0025\/2 with one roughness element (grid: 170 x 50 x 4). 
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Figure 86. Maximum amplitudes of the fundamental wave Ui, wj, mean-flow distortion uo, t>o, first har- 
monic wave u 2 , v 2 , and second harmonic wave u 3l v 3 for Re = 5000, k = 0.12 and e = 0.0025\/2 with two 
roughness elements (grid: 402 x 66 x 4). 
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perturbation streamfunction contours 




Figure 87. Contour plots of the instantaneous perturbation streaxnfunctions and vorticities for Re = 5000, 
k — 0.12, and e = 0.0025\/2 with seven roughness elements (grid: 402 x 66 x 4, flow direction is from left to right). 


total streamfunction contours 
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Figure 88. Contour plots of the total streamfunctions and vorticities for Re — 5000, k — 0.12, and e = 0.0025\/2 
with seven roughness elements (grid: 402 x 66 x 4, flow direction is from left to right). 
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Figure 89. M axi mum amplitudes of fundamental wave ttx, t>x, mean-flow distortion uo> ^o, first harmonic wave 
u 2 j and second harmonic wave U3, V3 for Re = 5000, k, = 0.12, and e — 0.0025\/2 with seven roughness 
elements (grid: 402 x 66 x 4). 
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9.2 Two-Dimensional Isolated Roughness Element for Flat Plate Boundary Layer 
Flow 

In this subsection, the effect of a two-dimensional isolated roughness element with different scales is studied. 
Since the reference scale is chosen to £q, the displacement thickness at inflow, the parameters for the surface grid 
and global grid generation also need to be changed. 

9.2.1 Grid Generation 

In this section, we use three kinds of roughness surfaces, which can be expressed as 

fi(x) ~ K/co3h 2 [y/b[{x - x 0 )], 

fl(x) = K,/cosh 2 [bi(x — Zo) 3 ]» 

/ 4 (z) = K/co$h 2 [bi(x - x 0 ) A ], (236) 

where k is the height of the roughness element, x 0 is the peak point x-coordinate, and bi is a parameter to adjust 
the length/height ratio. The grid in the physical (z,y) plane is not orthogonal, while in the computational (£,77) 
plane the grid is orthogonal and uniform. Noting that the computational domain in the normal direction is much 
large than that for channel flow, it is more reasonable to use a uniform grid in the physical plane above some 
height, yooi see Figure 90. The basic process of grid generation is: generate the grid in computational (£,77) plane; 
stretch it to an intermediate (xi,yi) plane; then curve it to fit the solid wall. 



(x,y) (*>Vi) ft' 7 *) 

Figure 90. Grid generation. 


The analytical mappings are 


V - 

yi 


{<? + ymax)yi 

<r + yi 

yoo*~±j> 0 < y < yoo, 


={ w 

l y 


3/00 y 1 /max* 


Thus, for y 6 [0, yoo], the resulting mapping is 


_ yoo(g+ Vma.)(y - /) 

^ 1/00(0' + y) — /(yoo + o') 


(237) 


or 

_ T?(y00 “ f) a + yOo/( cr + Vmax ~ V) 
y yoo(o- + Vmax ~ *l) 


( 238 ) 
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The required derivatives are: 


where 


Vx 




Vxx 




voo(g 4- Vmqg)(y - yoo )aj 1 
[yoo(ff + y) - /(yoo + <r)] 2 ’ 
yoo(g + yma«)(yoo - /)g 


yoo (»• + y) - /(yoo + o')] 2 ’ 
yoo(<^ + ymo*)(y - yoo)^ f 

/ _ . A #/ . M7 L, 


-SypoQ-Co- -t- y m «»)(yoo - /) 
[yoo(o- + y) - /(yoo + o')] 3 ’ 


(239) 


yoo 


(TVnj 0 

^ “I" 1/max tynjo 


(240) 


Usually, we choose njo = 2nj/3. For the y > yoo domain, the formulas are the same as for flat plate flow. Figure 
91 depicts a typical grid for this study. 



Figure 91. A typical grid used for simulation of rough flat plate boundary flow (part of the grid is shown) 


roughness size 

Re'o 

U) 


oli 

*0 

K 

small 

500 

0.070 

0.19418 

B222EHEHB 

35.41 

0.27 

medium 

663.12 

0.0968 

0.25970 

1 

44.0 

0.50 

large 

663.12 

0.0968 

0.25970 

-2.6244 x 10~ 3 

44.0 

1.12 


Table 4. Flow parameters for 2-D roughness flow. 
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roughness size 

C2d 

length 

T/At 

grids 

a 

Umax 

b 1 

small (fi) 

0.003 

(6+l)T-S 

400 

282 x 42 

4.5 

50 

0.25 

small (/a) 

0.003 

(6+l)T-S 

400 

282 x 42 

4.5 

50 

0.13 

medium (/1) 

0.004 

(6+l)T-S 

400 

226 x 42 

4.0 

50 

0.25 

medium (/a) 

0.004 

(6+l)T-S 

400 

226 x 42 

4.0 

50 

0.13 

large (/i) 

0.004 

(6+l)T-S 

400 

226 x 50 

4.5 

75 

0.25 

large (/<) 

0.004 

(6+l)T-S 

400 

226 x 50 

4.5 

75 

0.005 


Table 5. Computational parameters for 2-D roughness flow. 


9.2.2 Small-Scale Roughness 

We first consider a small roughness element with k = 0.27. The inflow section is chosen at Re* 0 = 500, the 
peak point of roughness is located at x 0 = 35.41 (measured from inflow section), and thus the local Re* = 550. 
The effective length of the roughness element defined as the distance between two / = 0.1k points, which is about 
3.64. 


2 /mai 
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Figure 92. Computational domain. 


Two different shapes of roughness, with the same roughness height k , are used to test the effect of both 
the size ratio and curvature of roughness elements on flat plate flow. In both cases, we also keep the effective 
length/height ratio 7 unchanged (see Tables 4 and 5). 

Case I: 

In this test, the shape of the roughness element is chosen to be 

fi(x) = 0.21 /cosh 2 [Q.b(x - 35.41)], 

where x is the nondimensionalized streamwise coordinate measured from the inflow boundary 
domain used here is 

x € [0,226.5], yE [0,50], 

which includes a six T-S wavelengths physical domain and a one T-S wavelength buffer domain. For this 2-D 
simulation, we use a 282 x 42 grid, which implies that the grid points for each streamwise T-S wavelength is 40. 

Due to the small dimensions of the roughness element and the small ratio of height/length, 7 (in this case, 
according to our definition, 7 = 0.074), no obvious separation bubble appears in the base flow, so, the flow recovers 


(241) 

. The computational 
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the Blasius profiles as soon as it passes the roughness sone. Figure 93 compares the computational results of the 
streamwise component of the base flow and the Blasius profile, showing no significant difference between these 
two profiles. 

It is nontrivial to describe the character of the local displacement 6*, momentum thickness 0, and shape factor 
H for this special flow, which are defined as 

• displacement thickness 


6 * 



- v)dy , 


• momentum thickness 


=r« 


(1 - |u|)dy, 


• shape factor 


H = 


6 ^ 
9 ' 


For Blasius flow, we have 


6* = 1.7208®/ v/feT, 0 = 0.664 x/y/Re^, (242) 

so, the shape factor is constant ( H « 2.59) throughout the plate. Figure 94 gives the streamwise distribution of 
£*, 9 and H } shows that H preserves the Blasius value except in the roughness region. 

A periodic disturbance with u> = 0.07 is then imposed at the inflow boundary to investigate the effect of 
roughness on disturbance amplification. The amplitude of the 2-D disturbance is set to 0.003, and each T-S 
period is divided into 400 time steps to ensure required time accuracy. The flow behaves similarly to that without 
roughness, but has a slightly enlarged amplification factor. Figure 95 depicts the behavior of the disturbance for 
both u and v. Except for a very small region, the amplitude of the disturbance grows smoothly. 

To be more precise, the inverse Fourier transformation is again employed here to investigate the spectral 
behavior of this flow. Compared to a smooth flat plate (see Figure 96), we find that the mean flow distortion, 
fundamental wave and first harmonic wave increased a little faster than the latter. Since the fundamental wave 
is two-orders larger in magnitude than the others, the major contribution of the enlargement comes from this 
component. The base flow and perturbation stream function contours are also given in Figure 97. 

Case II: 

Here we consider the roughness shape defined by 

f 2 (x) = 0. 2 7/cos/i 2 [0.13(a; - 35.41) 2 ]. 

All other parameters are the same, and the results are given in Figures 98-102. From Figure 98, we can see that 
the difference between the base flow and the Blasius profile is slightly enlarged as the shape of the roughness 
becomes more full. Also, from Figure 99, we see that the distortion region of the shape factor H is slightly larger 
than for Case I. Similarly, the amplification of the disturbance is slightly larger than for Case I, but there is little 
visual difference because of the similar base flow. 
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Figure 93. Comparison of the normal direction distribution of the streamwise velocity component of the 
base flow with small-scale roughness case I and the Blasius profiles at (a) x = 20, (b) x = 30, (c) x = 40, 
(d) x — 60. Roughness element located at x = 35.41, k = 0.27, ReJ = 500. 
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Figure 94. Streamwise distribution of and H for small-scale roughness case I. 
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Figure 95. Comparison of the amplification factor of disturbance for small-scale roughness case I and smooth flat 
plate. Re 5 = 500, k — 0.27, u = 0.07, e = 0.003. 
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Figure 96. Maximum streamwise and normal amplitudes of mean flow distortion uq, uo, fundamental wave tii, 
vi, and first harmonic wave ua, va, with nondimensional downstream distance for small-scale roughness case I 
and smooth flat plate. Re$ = 500, lj — 0.07, € — 0.003. 
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base flow streamfunction contours 
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Figure 97. Contour plots of steady base flow streamfunctions and instantaneous perturbation streamfunctions at 
t n 9 T for small-scale roughness case I. Flow direction is from left to right. 
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Figure 98. Comparison of the normal direction distribution of the streamwise velocity component of the base 
flow with small-scale roughness case II and the Blasius profiles at (a) x = 20, (5) x — 30, (c) x — 41, ( d ) i = 65. 
Roughness element located at x = 35.41, k = 0.27, ReJ = 500. 
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Figure 99. Streamwise distribution of 5*, 0, and H of the base flow for small-scale roughness case II. 




Figure 100. Comparison of the amplification factors of disturbance for small-scale roughness case II and smooth 
flat plate. Re 5 = 500, w = 0.07, e = 0,003. 
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Figure 101. Maximum streamwise and normal direction amplitudes for mean flow distortion tio, Vo, fundamental 
wave tii, vi, and first harmonic wave ua, va, with nondimensional downstream distance for small-scale roughness 
case II and smooth flat plate. J?eJ = 500, u ~ 0.07, e = 0.003. 
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9.2.3 Medium-Scale Roughness 


Case I: 

A slightly larger, single roughness element with k = 0.5 is considered here for the flat plate. The roughness 
element is located at x = 44.0 (peak point) and the inflow Reynolds number is Re* 0 = 663.12. The flow and 
computational parameters can be found in Tables 4 and 5, 

A separation bubble now clearly appears behind the roughness element (see Figure 103), with the length of 
the separation zone about 5-10 times the roughness height. The flow recovers the Blasius profile at about x = 80 
(thus the distortion region is about 60 times the roughness height ). Figure 104 depicts the comparison of the 
computed streamwise component of the base flow and the Blasius profile, showing that some significant changes 
occur near the roughness and that the flow almost recovers the Blasius profile at x = 85. Also, the streamwise 
distribution of the nondimensional displacement thickness 6 *(x), momentum thickness 0(z) and shape factor H(x) 
are given in Figure 105, showing that H is now clearly different from that of the Blasius flow in a much larger 
region than for the small roughness cases. 

An unstable mode is then imposed at the inflow boundary. The parameter is obtained from Dovgal & Kozlov 
(1990) so that some comparison can be performed later. In this case, according to our reference length and time, 
we have u> = 0.0968. The amplitude of the disturbance is set to e = 0.004, and one T-S period is still divided into 
400 time steps. 3600 time steps are used to obtain the final results. As shown in Figure 106, the amplification 
factor now is much larger than that for the smooth flat plate because of the base flow separation. 

By applying the inverse Fourier transformation, we decompose the flow field into the mean flow distortion wave 
u 0) uo, fundamental wave Ui, Vi, and first harmonic U 3 , v 3 ; see Figure 107. Obviously, all of these components 
are strongly amplified after passing the roughness element, with the high frequency wave showing the greatest 
increase. 

The instantaneous streamfunction contours for the perturbation flow at t = 9 T are given in Figure 108, 
showing the intensity of the disturbance increase. (Note that the final wavelength is in the buffer zone, so that is 
not correct.) 

Case II: 

Similar to the last section, we change the surface function to 

/ 2 (x) = 0.5/cos/i*[0.13(x — 44.0)*], (243) 

and keep all the other flow and computation parameters unchanged. The separation bubble is larger than that 
in Case I (see Figure 109), and the distortion region is also longer. From Figure 110, we find even at x = 85 that 
the flow has not fully recovered the Blasius profile. In this case, the length of the distortion region is about 80 
times the height of the roughness element. The distribution of 6*, 9 and H is also given in Figure 111. 

By introducing the same periodic disturbance as Case I, we obtain the amplification factor as depicted in 
Figure 112, which shows the amplification of the disturbance becoming stronger. Similar behavior can be found 
in Figure 113, in which the mean flow distortion, fundamental wave and first harmonic wave of the disturbance are 
amplified even more. Figure 114 gives the instantaneous perturbation streamfunction contour plots, also showing 
stronger amplification. 
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Figure 103. Streamfunction contours of the base flow for medium scale roughness case I. Flow direction is from 
left to right. 
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Figure 104. Comparison of normal direction distribution of the streamwise velocity component of the base flow 
with medium scale roughness case I and Blasius profiles at (a) x = 52, (b) x = 85. Roughness element location 
at x = 44.0, k = 0.5, Re* 0 = 663.12. 
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Figure 106. Comparison of the amplification factor of disturbance for medium-scale roughness case I 
Re* 0 = 663.12, w = 0.0968, e = 0.004. 
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Figure 107. Maximum streamwise and normal amplitudes of mean flow distortion uo, Uo> fundamental wave t*i, 
t/i, and first harmonic wave uj, V 2 } with nondimensional downstream distance for medium-scale roughness case I 
and smooth flat plate. Re$ = 663.12, u = 0.0968, e = 0.004. 
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Figure 108. Contour plots of instantaneous perturbation streamfunctions at t = 9T for medium-scale roughness 
case I. Flow direction is from left to right. 
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Figure 109. Contour plots of steady base flow streamfunction for medium-scale roughness case II. Flow direction 
is from left to right. 
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Figure 110. comparison of the normal direction distribution of the streamwise velocity component of the base 
flow with medium-scale roughness case I and the Blasius profiles at (a) x = 20, (6) x - 40, (c) x = 52, (d) x = 85. 
Roughness element located at x = 44.0, k = 0.50, ReJ = 663.12. 
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Figure 112. Comparison of the amplification factor of disturbance for medium-scale roughness case II and smooth 
flat plate. FeJ = 663.12, u = 0.0968, e = 0.004, 
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Figure 113. Maximum streamwise and normal amplitudes of mean flow distortion tio, i>o, fundamental wave 1 * 1 , 
Vi, and first harmonic wave U 2 , i> 2 , with nondimensional downstream distance for medium-scale roughness case 
II and smooth flat plate. = 663.12, w = 0.0968, e = 0.004. 
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Figure 114. Contour plots of instantaneous perturbation streamfunctions at t = 9 T for medium-scale roughness 
case II. Flow direction is from left to right. 


0.2.4 Large-Scale Roughness 

The experimental work of Klebanoff and Tidstrom (1972) and Dovgal fc Kozlov (1990) show that large-scale 
roughness elements located downstream will significantly change the amplification behavior of the disturbance. 
In this section, flow in the presence of a large isolated roughness element is studied. The parameters are given 
in Tables 4 and 5; for the second case, we use an even higher order wall function to approximate the rectangular 
roughness element used by Dovgal &; Kozlov (1990). 

Case I: 


The surface of the roughness is given by 

fi(x) = 1.12/cos/i a [0.5(x - 44.0)], (244) 

and other required parameters are given in Tables 4 and 5. 

The central difference scheme exhibited some trouble with this large scale roughness at the front of the 
hump. Instead of generating the first separation bubble, it only produce some oscillations, so the results are not 
physical. By applying the hybrid scheme developed by Patankar & Spalding (1972) to this problem, though both 
separation bubbles are generated, the rear separation region becomes much shorter than that produced by the 
central difference scheme. 


To obtain a reasonable base flow, a mixed scheme is then employed. From (224), we find that 

df 

— > 0 x < x 0) 
ox 


X > X 0 . 



118 


Liu^iUfMcCormick, Final Report 


Thus, the mixed scheme for the base flow is modifled to: 

hybrid: if ff > °, 

central difference: if < 0. 

Both separation bubbles in front and rear of the roughness element are then produced. The second bubble is 
about 20 times the height of the roughness element. Figure 115 gives the contour plots of the streamfunctions of 
the computed base flow. The comparison of the streamwise profile of our numerical results with the experimental 
results obtained by Dovgal & Koslov (1990) is depicted in Figure 116. Basically, they agree each other well for a 
large portion of the domain. The major difference appears upstream of the hump. The experiments show a larger 
separation zone ahead of the hump, but the computation gives a smaller separation zone, due apparently to the 
artificial viscosity introduced by the up-winding strategy and the surface shape difference between Case I and the 
experiment. The streamwise shape factor for our results and those of Dovgal & Koslov (1990) are also compared 
in Figure 117, which shows qualitatively good agreement. Note that all the geometric parameters appearing in 
the figures in this section have been interpreted according to those used in the experiments. 

A periodic perturbation with u> = 0.0968 (69Hz) and € =0.003 is then imposed at the inflow. In this case, one 
T-S period is divided to 400 time steps. Comparison of the computed amplification factors along the streamwise 
direction with those obtained by Dovgal & Kozlov (1990) for the hump is shown in Figure 118. Similar to the 
experiment, we obtained an exponential growth region downstream of the roughness element, which is mainly due 
to the effect of the rear separation bubble. Unfortunately, because the shape of roughness elements is different, 
the front separation bubble we obtained is smaller than that obtained by experiment for flow over a square hump. 
The disturbance almost has behavior similar to the linear instability stage before it reaches the rear separation 
bubble, and so, we do not fully obtain the first high rate growth upstream of the hump. The comparison of the 
amplitude of disturbance at select cross-sections is given in Figure 119, showing good qualitative agreement. 

Case II: 

Similar to the last subsection, we then changed the roughness shape so that it a better approximates the 
square roughness element. The shape function now is chosen to be 

f(x) = 1.12/ cosh 2 [0.005(as - 44.0) 4 ]. (245) 

Similar to Case 1, the mixed scheme is used to discretize the base flow. Both front and rear separation bubbles are 
larger than for Case I (see Figure 120). The shape-factor obtained is in quite good agreement with experimental 
results (Figure 121). The streamwise profiles also show better agreement with experiment (Figure 122). 

The amplification factors of the disturbance (Figure 123) and the normalized perturbation streamwise velocity 
profile (Figure 124) are also compared with those obtained by experiment, and both show better agreement 
than for case I. The perturbation streamfunction contours are given in Figure 125. Note that some additional 
interactions between the experimental measurement instrument and the base disturbance level of the wind tunnel 
may generate some 3-D disturbance, which is not possible for this 2-D simulation to describe. 
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Figure 115. Streamfunction contours of the base flow for large-scale roughness case I. Flow direction is from left 
to right. 



Figure 116. Comparison of the numerical streamwise velocity profile for large-scale roughness case I and experi- 
mental results obtained by Dovgal & Kozlov (1990, over a square hump). 
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Figure 117. Comparison of the numerical streamwise disturbance of shape factor H and the experimental results 
obtained by Dovgal & Kozlov (1990, over a square hump). 
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Figure 118. Streamwise distribution of the amplification factors: (a) numerical results for large-scale 
roughness element, (b) experimental results obtained by Dovgal & Kozlov (1990, over a square hump). 
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Figure 119. Normalized amplitude of streamwise perturbation velocity: (a) numerical results for 
large-scale roughness case I, (b) experimental results obtained by Dovgal & Kozlov (1990, over a square hump). 
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Figure 120. Contour plots of the base flow for large-scale roughness case II. 
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Figure 121. Comparison of numerical shape factor distribution for large-scale roughness case II and experimental 
results obtained by Dovgal & Kozlov (1990) for a square hump. 
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Figure 122. Comparison of streamwise velocity profile at selected locations with experimental results. 
Re* 0 = 663.12, k = 1.12. 
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Figure 123. Comparison of amplification factors: (a) large-scale roughness case II, 
(b) experimental results obtained by Dovgal & Kozlov (1990) over a square hump. 
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Figure 124. Comparison of the normalized streamwise disturbance velocity profiles: (a) numerical results 
for large-scale roughness case II, (b) experimental results of Dovgal & Kozlov (1990) over a square hump. 
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Figure 125. Contours plots of instantaneous perturbation streamfunctions at t = 9 T, 
Jtej = 663.12, u = 0.0968, e = 0.003. Flow direction is from left to right. 


9,3 Two-Dimensional Distributed Multiple Roughness 

Since the major contribution to the disturbance amplification comes from the separation bubbles, it is expected 
that the multiple roughness will induce larger amplification factors through the whole domain. To test this, five 
roughness elements are now chosen, described by 

f(x) = 0.5/co$/i 2 [0.005(x - 44) 4 ] + 0.7 /cosh 2 [0.005(b - 79) 4 ] + 0.9/cos/i 2 [0.005(x - 114) 4 ] 

+ 0.7/cos/i 2 [0.005(x - 149) 4 ] + 0.5/cosh 2 [0.005(x - 184) 4 ]. (246) 

The computational domain is 


xe [0,242], ye [0,50], (247) 

the grid stretch parameter o = 4.5, and the grid sise is 362 x 50 (9 T-S wavelengths physical domain and 1 
T-S wavelength buffer). The contour plots of the base flow are depicted in Figure 126, and the associated shape 
parameter H and both the displacement and momentum thickness are shown in Figure 127. 

For the perturbation flow, the parameters are the same as for the large roughness cases, but the amplitude 
of disturbance is reduced to e=0.0005. As expected, the amplification factors of the disturbance reach a higher 
value. Since the roughness elements are distributed quite randomly, the fundamental T-S wave is almost invisible 
after passing a few roughness elements. Figure 128 displays the amplification factors for both u and if, and Figure 
129 gives the instantaneous streamfunction contours at t = 10T. 



Figure 127. Nondimensional shape parameter H , displacement thickness 6* and momentum thickness 9 of the 
base flow with the presence of multiple roughness elements. 
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Figure 128. Streamwise distribution of amplification factors; 
(a) amplification factors for u, (b) amplification factors for v . 
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Figure 129. Contour plots of the instantaneous perturbation streamfunctions at t = 10 T for the distributed 
multiple roughness elements (part of the domain is shown, flow direction is from left to right). 
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9.4 Three-Dimensional Isolated Roughness Element for Flat Plate Boundary Layer 
Flow 


In this subsection, the effect of a three-dimensional isolated roughness element with different scales is studied. 
The parameters for this study are listed in Tables 6 and 7. 


The computational domain is set to 

X G [0, mL ,], y G [0, y m „], z G [-^, |], 

and the shape of the surfaces of the roughness is defined by 

f{x, z) = K/cosh 2 [bi(x - X 0 ) 4 + M 4 )]. 

The roughness Reynolds number 



v 


where U K is the undisturbed velocity at height k, are also listed in Thble 6. 


(248) 


size 

Re'o 

Re* 

U) 


«/ 

(3 

*0 

K 

small 

663.12 

31.83 

0.0968 

0.25970 

-2.6244 x NT 3 

0.25 

44.0 

0.30 

medium 

663.12 

123.34 

0.0968 

0.25970 

-2.6244 x 10- 3 

0.25 

44.0 

0.60 

large 

663.12 

418.14 

0.0968 

0.25970 

-2.6244 x 10- 3 

0.25 

44.0 

1.12 


Table 6. Flow parameters for 3-D roughness flow. 


size 

*2d 

length 

T/At 

grids 

<7 

Vmax 

bi 

b 2 

small 

0.003 

(3+l)T-S 

360 

98 x 34 x 34 

4.0 

50 

0.005 

0.0025 

medium 

0.003 

(3+l)T-S 

360 

98 x 34 x 34 

4.0 

50 

0.005 

0.0025 

large 

0.003 

(3+l)T-S 

400 

114 x 34 x 34 

4.0 

50 

0.005 

0.0025 


Table 7. Computational parameters for 3-D roughness flow. 


0. 4.1 Base Flow 

The base flow is calculated by the hybrid scheme, which yields better stability and faster convergence than 
central differences, though less accurate for large mesh Reynolds number. Three cases are under consideration. 

1. Small-Scale Roughness In this case, we set k = 0.3 and the roughness Reynolds number Re K = 31.83. 
Other parameters can be found in Tables 6 and 7. Though the base flow is still streamwise dominant, it is really 
three-dimensional, i.e., the spanwise velocity component, w, is no longer sero. Figure 130 depicts the velocity 
vector plots of the base flow field on the j = 2 and j = 3 grid surfaces. No obvious separation is observed at 
either the front or rear of the roughness element, because the roughness surface function that we used is quite 
smooth and the roughness height is really small. 

2. Medium-Scale Roughness We then change the roughness height to k = 0.6, and the corresponding 

roughness Reynolds number to Re* = 123.3 < 300. A separation sone is found behind the roughness, which is 
much smaller than what we observed in the 2-D roughness case with the same roughness Reynolds number. The 

vector plots of the base flow on j = 2 grid surface is show in Figure 131, demonstrating that a very small separation 
sone exists in front of the roughness and another separation sone appears at the rear of the roughness. The front 
separation can form a horseshoe-shaped filament, which wraps around the front of the roughness element, as it 
travels downstream and is further increased. 
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3. Large-Scale Roughness As the height of the roughness /c is increased to 1.12 (the same as that used for 
2-D case), the roughness Reynolds number increases to Re K = 418.14 > 300. Two separation bubbles are clearly 
found in front and at the rear of the roughness element. The front separation forms a horseshoe-shaped vortex 
filament (see Figure 132 for the contour plots of streamwise vorticity on selected (y, z)- planes), and may cause 
a secondary instability when it travels downstream. Besides, a counter-rotating pair of vertical spirals is found 
close to the rear wall and near the top of the roughness element. Figure 133 gives the vector plots of the velocity 
field on the j = 3 grid surface. The magnitude of vorticity contours on the y = 0.245 plane are depicted in Figure 
134, which clearly exhibits the existence of the pair of vertical spirals. 

This base flow pattern obtained by our computation qualitatively supports the experimental results obtained by 
Mochizuki (1961a,b). However, Mochizuki used a sphere as the isolated roughness element, and assumed elements 
of nearly equal height, span, and length, which is very different from ours. In our examples, the length and span 
are several times larger than the height. Therefore, it is difficult to obtain any quantitative comparison between 
our computational results and Mochizuki’s experiment; yet the qualitative agreement between computation and 
experiment is nonetheless encouraging. 


9.4.2 Disturbance 

After obtaining the base flows numerically, we then impose a 2-D periodic small disturbance (T-S wave) at the 
inflow boundary and study the behavior of the developing perturbation. Again, three cases described in Tables 
6 and 7 in the last subsection are investigated. 

1 . Small-Scale Roughness Figure 135 depicts the magnitude of disturbance vorticity contours on the j = 2 
and j = 3 grid surfaces, which clearly shows that the disturbance has become three-dimensional, even though no 
3 -D disturbance is imposed at inflow. This is different from the case of 2-D roughness elements. Figure 136 gives 
the amplification curve for the velocity magnitude of perturbation and its u— and w— components. Since the 
flow has become three-dimensional, we use the amplification factor of disturbance velocity magnitude given by 

maz{v/ii /2 + v ' 2 + u/ 2 } 

Ny — ln( 7 = — ) 

max{ y + u#} 

as the measure of amplification. Here, U * , v * and w * are perturbation velocity components, and the subscript 0 
refers to inflow values. A similar method can be used for the 2-D problem. From Figure 136 we find that when the 
flow meets the roughness element, the flow perturbation gains more energy, but past the roughness zone energy 
decreases. 

2. Medium-Scale Roughness Evidently, In this case, the disturbance is enlarged much faster than for the 
small-scale roughness case, and the disturbance in the spanwise direction grows very rapidly. Figure 137 gives the 
contour plots of the magnitude of disturbance vorticity on the j = 2 and j — 3 grid surfaces. The fundamental 
2 -D T-S wave is apparently stimulated by streamwise vorticity carried by the travelling horseshoe vortex. 

3 . Large-Scale Roughness The original version of the code has difficulty with disturbances in the case of 
large-scale roughness elements. Because of the dramatic amplification of the disturbance, central differences 
induce the flow to break up in front of the roughness element. Although the separation area is much smaller than 
that for 2-D roughness, the flow more easily breaks down for 3-D roughness. It seems that the formation and 
development of the horseshoe vortex strongly undermines flow stability. In this computation, though the major 
disturbance has not passed the roughness region, we still observe some 3-D disturbance behind the roughness 
element. Thus, it seems that the behavior of the disturbance is changed not only in space, but also in time. Some 
disturbance with greater travelling speed is generated. 
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Figure 130. Vector plots of the base flow on j = 2 grid surface for the small-scale 
3-D roughness case. Re* 0 = 663.12, Re K = 31.83, k = 0.3. 



Figure 131. Vector plots of the the base flow on j = 2 grid surface for the medium-scale 
3-D roughness case. i2ej = 663.12, Re K — 123.3, /c = 0.6. 



Figure 132. Contour plots of the streamwise vorticity (base flow) on i = 45, 50, 55, 60, 65, 70,75, 80 
(y, z)-planes for the large-scale 3-D roughness case. i2ej = 663.12, Re K = 418.14, k = 1.12. 
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Figure 133. Vector plots of the base flow on j = 3 grid surface for the large-scale 
3-D roughness case. = 663.12, Re K — 418.14, k = 1.12. 



Figure 134. Contour plots of the spanwise vorticity (base flow) on k = 12, 17, 22 (x,y)-planes 
for the large-scale 3-D roughness case. = 663.12, Re K — 418.14, k = 1.12. 



u 0 



Figure 135. Contour plots of the disturbance vorticity magnitude on j — 2 and j = 3 grid surface (from 
bottom to top) for the small-scale 3-D roughness case at t = 5 T. = 663.12, k — 0.3, e = 0.003. 
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Figure 136. (a) Comparison of amplification factors of u between the 2-D smooth fiat plate and 3-D small 
-scale isolated roughness case, (b) Comparison of the disturbance velocity magnitude amplification factors 
between 2-D smooth flat plate and 3-D small-scale isolated roughness case, (c) Amplification factors of 
spanwise velocity magnitude w for the 3-D small roughness case. Rcq = 663.12, k = 0.3, c = 0.003. 
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Figure 137. Contour plots of the perturbation vorticity on (a) j — 2 and (b) j = 3 grid 
surfaces, for the medium-scale 3-D roughness case. — 663.12, k = 0.6, e = 0.003. 


9.5 Three-Dimensional Distributed Roughness Element for Flat Plate Boundary 
Layer Flow 

The multiple roughness case is as a collection of single three-dimensional roughness elements distributed 
randomly in the domain, with arbitrary height and spacing. According to Morkovin and Corke (1984), distributed 
roughness will cause 

• complex, predominantly streamwise vorticity patterns; 

• long induction distances from the process-initiating roughness to the onset of turbulence. 


General multiple roughness simulation is not yet possible. In this section, a simplified model with seven 
roughness elements in the computational domain is used (see Figure 138). By assuming the spanwise wavenumber 
(3 — clr, the computational domain is then restricted to 

x £ [0, 121], y £ [0, 50], z £ [-12.1, 12.1], 

and the shape of the surfaces with seven roughness elements is defined by 

f(x,z) = K/cosh 2 [b{x — 29. 8) 4 4* bz 4 )] + K/cosh 2 [b(x — 54. 0) 4 + bz 4 )] + 

k/ cosh 2 [b(x — 78. 2) 4 -f bz 4 )] + K,/cosh 2 [b(x — 41. 9) 4 -f b(z — 12. 1) 4 )] + 

K/cosh 2 [b(x — 41. 9) 4 + b(z *f 12. l) 4 )] + K,/cosh 2 [b(x — 66. 1) 4 + b(z — 12. l) 4 )] + 
K,/cosh 2 [b(x — 66. 1) 4 + b(z -f 12. l) 4 )]. 

A 122 x 34 x 34 grid is used, which includes a 4.3 wavelengths physical domain and a 0.7 wavelength buffer 
domain. 
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9.5.1 Base Flow 

Two cases with different sizes of roughness are considered. 

1. Small-Scale Distributed Roughness 

First, let b = 0.005 and k — 0.35 [Re K ~ 43.5), which is considered to be small-scale roughness. Figure 139(a) 
depicts the vector plots of velocity field on the j = 2 and j — 3 grid surfaces). Yet the base flow is no doubt three- 
dimensional with a w— component. In addition, the so-called horseshoe vortices are generated downstream of the 
computational domain, which is different from the isolated roughness flow. Figure 139(b) gives the streamwise 
vorticity contours on the selected (y, z)-planes, clearly shows the development of horseshoe vortices. 

2. Large-Scale Distributed Roughness 

We then increase the height of roughness to k — 1.12, yielding a corresponding roughness Reynolds number is 
of Re K ~ 418, which is considered to be large-scale roughness. Similar to what we described in the case of a single 
3-D roughness element, two separation zones in front and at the rear of each roughness element are found. The 
rear are stronger than the front vortices. The rear vortex energy comes from the downwash of the flow from the 
high velocity fluid at the top of the roughness element (high momentum region), while the front horseshoe vortex 
energy comes from stagnation of the fluid in front of the roughness element (low momentum region). Figure 
140(a) depicts the vector plots of velocity field on the j = 2 and j = 3 grid surfaces, and Figure 140(b) gives 
the streamwise vortex contours on the selected (y, z)-planes, clearly showing the strong development of horseshoe 
vortices. 

Because of the separation at the rear of the roughness elements, the profiles of the base flow are distorted from 
the Blasius solution. Figure 141 depicts the u— velocity profiles between three roughness elements (in the central 
(z, y)-plane of the computational domain), which shows qualitatively agreement with experimental measurements 
by Tadjfar et al.(1993). 


0.5.2 Disturbance 

We then impose a 2-D periodic perturbation (T-S wave) on the base flow. Unfortunately, the central difference 
scheme shows numerical instability for the large-scale roughness case with early breakup. On the other hand, 
the code works well for the small-scale roughness case. With = 0.005, and one T-S period divided into 360 
time steps, we obtained the expected results. Figures 142 and 143 depict the amplification curve of the maximum 
perturbation velocity and its u— component, which is several times larger than what was obtained for the 2-D 
T-S wave past the smooth flat plate around the roughness region. According to the experiments of Leventhal 
& Reshotko (1981), for flow over multiple roughness elements, the maximum amplification rates are about three 
times as large as those for T-S instability over a smooth plate. Our computation qualitatively supports their 
experiments. Figure 144 gives the contours of spanwise vorticity on the j — 2 and j = 3 grid surfaces, which 
shows that the T-S waves are apparently weakened when the disturbance passes the distributed roughness zone, 
while the A— wave gradually develops downstream, and thus the flow stability worsens. This is also qualitatively in 
agreement with the experiments of Tadjfar et al. (1993) in which the oscillation of disturbance amplified rapidly 
as it moved downstream, although no evidence of T-S waves in the boundary layer was observed. 

To investigate the effect of the amplitude of the imposed disturbance, we then increased the amplitude of 
disturbance € 2 d — 0.015. The H-type like pattern now can be seen clearly (see Figure 145). When the intensity 
of the disturbance is increased further (e 2 d = 0.05), the hairpin vortex (cf. Morkovin, 1990) clearly appeared (see 
Figure 146), and the laminar flow then broke down. 
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Figure 138. Distribution of the 3-D distributed roughness surface. 



Figure 139(a). Vector plots of velocity fields on j = 2 and j — 3 
grid surfaces for the small-scale 3-D distributed roughness case. 
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Figure 140(b). Contour plots of streamwise vorticity on i = 30,42, 54, 66, 78, 90 
(y, <z)-planes for the large-scale 3-D distributed roughness case. 



(a) 



(b) 

Figure 141. Comparison of u— velocity profiles for the base flow, (a) Numerical results for the large-scale 
3-D distributed roughness case, (b) experiment of Tadjfar et al (1993, over the sphere roughness). 
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Figure 142. Amplification factor of the maximum perturbation velocity for the small-scale 3-D distributed 
roughness obtained on a 82 x 34 x 18 (4.4 T-S wavelengths physical domain + 0.6 wavelength buffer) grid. 
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Figure 143. Amplification factor of the perturbation u-velocity for the small-scale 3-D distributed roughness 
obtained on a 82 x 34 x 18 (4.4 T-S wavelengths physical domain -b 0.6 wavelength buffer) grid. 


Figure 144. Contours of spanwise perturbation vorticity on j = 3, 5 and 7 
grid surfaces (from bottom to top) for the small-scale distributed roughness. 
e 2 d — 0.005, k = 0.35, Re K = 43.5, flow direction is from left to right. 



Figure 145. Contour plots of relative helicity of the total flow at 
t — 5T for the small-scale distributed roughness. e 2 d — 0.015, 
k — 0.35, Re K = 43.5, flow direction is from left to right. 




Figure 146. Contour plots of relative helicity of the total flow at 
t = 5T for the small-scale distributed roughness. € 2 d = 0.05, 
k = 0.35, Re K — 43.5, flow direction is from left to right. 
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10 Concluding Remarks 

A very accurate and efficient multigrid code has been developed in Phase II of the project. This new code 
supports direct numerical simulation for 3-D smooth channels, flat plates, and flows with boundary layers incor- 
porating 2-D and 3-D single and multiple roughness elements. 

• The fourth-order and fully-implicit time-marching scheme on a staggered and stretched grid is highly accu- 
rate and stable for transition problems. 

• The semi-coarsening method based on line-distributive relaxation is very efficient, making the cost of the 
fully-implicit scheme comparable to explicit schemes at each time step. 

• The new version of the governing equations in general coordinates, which considers all u , v , tu, P, 17, V, W 
as unknowns, together with a high-order mapping successfully performs DNS for transitional flow with 
roughness elements, and is extendible to problems with more general geometry. 

• The new treatment of the outflow boundary dramatically reduces computational cost, making realistic 
spatial simulation much more feasible. 

• The surface roughness contributes to early transitions in different ways for 2-D or 3-D single or multiple 
roughness elements. The current DNS computation shows quantitative agreement with experiments by 
Dovgal and Kozlov (1990) for a single 2-D roughness elements and qualitatively agreement with experiments 
by Mochizuki (1961a,b) for a single 3-D roughness element and by Tadjfar et al. (1993) for multiple 3-D 
roughness elements. 
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